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Nonlinear,  Dispersive  Long  Waves  in 
Water  of  Variable  Depth 

James  T.  Kirby 

Center  for  Applied  Coastal  Research 

University  of  Delaware,  Newark,  DE  19716,  U.S.A. 

Abstract 

The  use  of  weakly-dispersive  models  to  compute  the  propagation  of  nonlinear 
ocean  surface  waves  in  the  coastal  environment  is  reviewed.  Models  which  are 
fully  two-dimensional  in  horizontal  coordinates  are  discussed  first,  and  vari¬ 
ous  approximations  involved  in  various  models  are  developed  and  compared. 
Available  tests  of  the  accuracy  of  model  predictions  are  also  reviewed.  At¬ 
tempts  to  incorporate  realistic  effects  including  wave  breaking,  shoreline  runup 
and  wave-current  interaction  into  model  schemes  are  discussed.  Subsequently, 
models  for  one-directional  and  weakly-two-directional  propagation  are  briefly 
reviewed.  Finally,  frequency-domain  formulations  are  reviewed.  In  all  cases, 
the  primary  concentration  is  on  aspects  of  model  development  related  to  the 
computation  of  realistic  waves  in  the  ocean  environment. 


1  Introduction 

Water  waves  in  shallow  water  (characterized  by  depths  which  are  small 
compared  to  a  wavelength)  deviate  more  severely  from  the  basic  assump¬ 
tions  of  linear  wave  theory  than  do  waves  in  deeper  water.  Waves  in 
shallow  water  are  modified  rapidly  over  distances  comparable  to  a  wave¬ 
length,  due  to  weak  frequency  dispersion  and,  as  a  result,  nearly  resonant 
interaction  between  sets  of  three  Fourier  harmonics.  These  resonant  in¬ 
teractions  lead  to  a  strong  transfer  of  energy  from  fundamental  to  higher 
harmonic  components  as  waves  approach  shore  over  natural  beaches, 
leading  to  waves  having  relatively  narrow  crests  and  broad  troughs.  In 


addition,  the  rapid  transfer  of  energy  from  low  to  high  frequencies  is 
manifested  in  front-to-back  asymmetry  of  the  wave  profile  as  the  wave 
pitches  towards  the  shore.  Finally,  difference  interactions  lead  to  a  rapid 
growth  in  low-frequency  energy  as  well,  leading  to  shore-normal  surf  beat 
or  longshore  propagating  edge  wave  and  leaky  mode  components. 

Linear  theory  does  not  give  us  much  of  a  route  to  describing  the  com¬ 
plexity  of  wave  motion  in  the  final  stages  of  shoaling  or  in  the  surf  zone. 
For  this  reason,  a  body  of  theory  and  computational  techniques  has  de¬ 
veloped  over  the  years  that  is  specifically  aimed  at  describing  the  regime 
where  frequency  dispersion  effects  are  small  and  nonlinearity  takes  on 
an  equal  or  more  important  role.  The  first  steps  in  the  study  of  these 
phenomena  were  made  by  Airy  [1],  who  obtained  the  Nonlinear  Shallow 
Water  (NSW)  equations  governing  nondispersive  nonlinear  wave  propa¬ 
gation,  and  by  Boussinesq  [2]  and  Korteweg  and  deVries  [3],  who  provided 
models  incorporating  the  leading  order  effect  of  frequency  dispersion  in 
the  long  wave  approximation.  This  body  of  work  has  played  a  central 
role  in  the  evolution  of  the  mathematical  physics  of  waves  that  exhibit 
a  balance  between  nonlinear  and  dispersion  effects.  In  particular,  the 
Korteweg-deVries  (KdV)  equation  has  been  proposed  as  the  principle 
model  equation  for  a  wide  range  of  physical  phenomena,  and  has  gener¬ 
ated  a  vast  body  of  literature  related  to  its  exact  solution  for  arbitrary 
initial  conditions  by  means  of  the  Inverse  Scattering  Transform  [4]  [5]. 
Many  textbook  treatments  of  this  basic  aspect  of  the  theory  for  waves  in 
an  infinite  or  periodic  domain  now  exist  [6]  [7]. 

In  the  open  coastal  environment,  the  physical  distance  traversed  by 
a  surface  wind  wave  propagating  from  the  point  where  it  begins  to  be 
significantly  affected  by  the  bottom  until  it  reaches  shore  is  often  no  more 
than  several  wavelengths.  Thus,  the  problem  of  computing  the  evolution 
of  waves  through  domains  with  significant  inhomogeneities  (occuring  over 
length  scales  which  are  not  terribly  long  compared  to  a  wavelength)  is 
of  central  importance.  The  problem  of  rapid  wave  evolution  in  a  varying 
domain  is  different  from  the  more  academic  problem  of  evolution  over 
long  distance  in  a  uniform  domain,  and  remains  more  a  problem  for 
numerical  analysis.  It  is  this  aspect  of  the  theory  for  weakly-dispersive 
long  waves  that  is  primarily  reviewed  here. 

The  modern  birth  of  the  study  of  nonlinear  wave  propagation  over 
an  uneven  bottom  took  place  in  the  late  1960’s  [8]  [9]  [10].  Peregrine  [9] 
provided  the  derivation  of  a  Boussinesq-type  model  using  depth-averaged 


velocity  as  the  dependent  variable,  and  provided  the  first  numerical 
calculations  of  the  phenomenon  of  nonlinear  wave  shoaling  and  evolu¬ 
tion.  Development  of  these  techniques  continued  through  the  1970’s  and 
80’s  [11]  [12]  [13],  but  awaited  the  development  of  both  theoretical  tech¬ 
niques  for  extending  the  range  of  applicability  of  the  models  (through 
improved  linear  dispersion  properties)  as  well  as  the  development  and 
availability  of  computers  capable  of  making  model  computations  on  a 
more  routine  basis.  Both  of  these  hurdles  have  been  crossed,  and  there 
has  been  a  resulting  explosion  of  interest  in  shallow  water  wave  models 
in  the  past  ten  years. 

In  this  review,  we  first  concentrate  on  the  development  of  the  Boussi- 
nesq  model  framework,  and  review  recent  advances  which  have  extended 
the  applicability  of  these  models  over  much  of  the  range  of  depths  where 
waves  feel  the  bottom.  Efforts  to  extend  the  Boussinesq  model  approach 
to  cover  wave  breaking,  shoreline  runup  and  wave-current  interaction  are 
also  reviewed.  Recent  results  in  the  area  of  reduced- dimension,  forward- 
propagation  models  are  also  considered,  and  we  finish  with  a  discusson 
of  spectral  approaches  to  wave  computations. 


2  Governing  equations  and  scaling  funda¬ 
mentals 

We  concentrate  here  on  a  review  of  the  development  of  models  for  the 
propagation  of  waves  through  a  homogeneous  inviscid  fluid  with  density 
p.  We  employ  a  Cartesian  coordinate  system  with  2  oriented  upward  and 
with  the  x,  y  plane  lying  in  the  plane  of  the  still  water  free  surface.  The 
fluid  is  bounded  below  by  a  surface  Ft,  =  z  +  h(x,  y,  t)  =  0  and  above 
by  a  surface  Fs  =  z  —  r)(x,y,t)  =  0.  Time  dependence  in  the  surface  is 
an  obvious  manifestation  of  wave  motion.  We  allow  for  time  dependence 
in  the  bottom  bounding  surface  to  allow  for  the  possibility  of  tsunami 
genesis  through  earthquake  motions  and  other  shifts  in  bottom  geometry. 

In  section  2.4,  we  will  introduce  scaling  assumptions  that  make  a 
strong  distinction  between  the  relative  size  of  horizontal  and  vertical 
motions.  For  now,  we  retain  a  more  general  formulation.  Introducing  a 
velocity  vector  u  given  by 

U  =  (u(x,Z,*),»(x,2,t),tfl(x,Z,t))  (1) 
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where  x  =  ( x,y )  is  the  position  vector  in  the  horizontal  propagation 
space,  we  may  write  the  governing  equations  as 

V3  •  u  =  0  (2) 

ut  +  u  •  V3«  +  ^V3p  +  giz  =  0  (3) 

where  V3  denotes  a  gradient  in  three  dimensions.  We  will  often  make 
a  distinction  between  horizontal  and  vertical  velocity  components.  We 
define  the  horizontal  velocity 

u  =  (u(x,z,t),v(x.,z,t))  (4) 

and  denote  a  gradient  in  horizontal  coordinates  by  V.  Kinematic  con¬ 
straints  on  the  bounding  upper  and  lower  surfaces  lead  to  the  conditions 

T)t  +  u  •  V?7  =  w,  z  =  ?/(x,  t )  (5) 

and 

ht  +  u  •  V/i  —  —  w,  z  =  —  h(x,  t)  (6) 

Finally,  we  ignore  dynamic  interaction  with  an  overlying  air  layer  and 
simply  specify  that  pressure  be  constant  on  the  upper  bounding  surface, 

P  =  0;  z  =  T](x,t)  (7) 

2.1  Irrotational  flow 

The  formulation  above  is  completely  adequate  for  the  general  problem 
of  inviscid,  incompressible  flow.  Several  of  the  formulations  considered 
below  will  rest  only  on  the  formulation  above,  without  recourse  to  further 
assumptions.  However,  due  to  the  implications  of  Kelvin’s  theorem  for 
inviscid,  barotropic  fluids,  waves  are  often  assumed  to  fall  into  the  cat¬ 
egory  of  irrotational  flows;  i.e.,  flows  having  no  vorticity  or  circulation. 
In  this  case,  stronger  conditions  may  be  imposed  on  the  formulation  of 
the  wave  motion  problem.  For  the  case  of  irrotational  motion, 

V3  x  u  =  0  (8) 

we  may  introduce  a  velocity  potential  <^(x,  z,  t)  such  that 

u  —  V3<£  (9) 
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Subsequent  use  of  (9)  in  (2)  then  gives 

Vl*=  0  (10) 

which  is  Laplace’s  equation  for  the  velocity  potential.  This  result  is 
supplemented  by  the  first  integral  of  the  Euler  equation  (3),  known  as 
Bernoulli’s  equation, 

-  +  <t>t  +  gz  +  -|V3<^|2  =  C(t)  (11) 

P  2 

which  relates  the  pressure  field  to  the  velocity  field.  Evaluating  (11)  at 
the  free  surface  and  introducing  the  constraint  (7)  leads  to  a  dynamic 
free  surface  boundary  condition 

9V  +  <f>t  +  ^|V3<£|2  =  0  (12) 

where  the  Bernoulli  constant  is  absorbed  under  the  assumption  that  the 
free  surface  conincides  with  the  level  surface  z  =  0  when  no  wave  activity 
or  current  field  is  present.  Finally,  the  kinematic  boundary  conditions 
are  given  by 


r)t  +  V<£  •  V»7  =  <f>z  ;  z  =  r 7 

V  (f> -Vh  — —<j)z  ;  z  =  —h 

2.2  Variational  formulation 


(13) 

(14) 


The  existence  of  a  Lagrangian  for  irrotational  free  surface  flows  was  in¬ 
troduced  originally  by  Luke  [14],  who  showed  that  the  Lagrangian  was 
given  by  the  integral  of  fluid  pressure  over  the  depth;  i.e., 

£(x,t)  =  f  p(-x,z,t)dz  (15) 

J —h 

Using  Bernoulli’s  equation  (11)  above,  we  may  recast  C  in  the  form 


£(x>  t)  =  \pgh2  -  P  f_h  <!>tdz  -  ^pgri2  +  P-  /j  W)2  +  {(f>z)2]dz  (16) 


Luke’s  principle  then  states  that  the  integral  of  C  over  all  horizontal 
space  and  time  is  then  stationary  with  respect  to  variations  in  t]  and  (j) 

x,  i;77,  <f>)dxdt  =  0 


(17) 


where  we  neglect  the  influence  of  lateral  boundaries  on  the  specification 
of  the  problem.  It  may  be  shown  [14]  that  this  formulation  is  equivalent 
to  the  boundary  value  problem  specified  by  equations  (10  -  14).  An 
examination  of  (16)  shows  that  the  first  term  is  trivial,  as  it  does  not 
depend  on  a  dynamic  variable.  The  last  term  in  brackets  consists  of  the 
sum  of  the  potential  energy  per  unit  area  in  the  horizontal  plane  V  and 
the  kinetic  energy  per  unit  area  in  the  horizontal  plane  T.  We  denote  the 
sum  of  these  two  energies  as  the  total  energy  W(x,t).  The  Lagrangian 
(16)  is  fully  equivalent  to  the  usual  Lagrangian  from  classical  mechanics, 

C  =  T  —  V  (18) 

This  result  may  be  established  by  use  of  Green’s  first  identity  and  Laplace’s 
equation  applied  to  the  kinetic  energy  integral  in  (16).  Following  Miles  [15], 
we  may  integrate  the  second  term  in  (16)  by  parts  once  and  obtain  the 
relation 

I-h^tdz  =  §i  J-h^dZ  ~  ^X’^;  £CM)  =  ^(x,z  =  »/,<)  (19) 

where  the  first  term  integrates  to  the  boundaries  in  time  in  (17)  and 
hence  has  no  dynamic  significance.  The  Lagrangian  may  then  be  written 
as 

£  =  p(nt  -  H  (20) 

This  result  may  be  further  identified  as  the  Legendre  transform  [16],  with 
~H  denoting  the  Hamiltonian,  the  canonical  coordinate,  and  p(  the  con¬ 
jugate  momentum.  Proof  that  the  surface  wave  problem  for  irrotational 
motion  is  a  Hamiltonian  system  was  provided  originally  by  Zakharov  [17]. 


2.3  Computational  approaches  using  the  full  equa¬ 
tions 

Given  the  full  boundary  value  problem  for  fluid  motion,  it  is  possible 
to  construct  numerical  solutions  for  the  full  three-dimensional  problem. 
This  approach  has  been  pursued  by  a  number  of  authors,  principally 
using  a  technique  known  as  the  Volume  of  Fluid  (VOF  method).  This 
method  has  been  applied  to  both  shoaling  waves  and  breaking  waves 
with  some  success  [18].  Alternately,  for  the  irrotational  problem,  several 
approaches  exist  which  reduce  the  dimensionality  of  the  problem  being 
studied  by  one.  Two  such  approaches  are  reviewed  here. 


2.3.1  Hamiltonian  canonical  evolution  equations 


The  existence  of  a  Hamiltonian  for  an  irrotational  flow  with  a  free  surface 
provides  an  immediate  path  to  developing  evolution  equations  for  wave 
motion,  in  the  form  of  the  canonical  equations  for  the  evolution  in  time 
of  the  coordinate  and  conjugate  momentum, 

it  =  (21) 

(.  =  ~{%()  (22) 

where  the  derivatives  are  functional  derivatives  with  respect  to  the  de¬ 
pendent  variables.  Following  Miles  [15],  (21)  and  (22)  may  be  written 
as 


Vt  =  -Ve-Vif  +  CU  +  ViyVif)  (23) 

(,  =  -S-)-iv$-V{+ic2(l  +  V,.V,)  (24) 

where  (  =  </>2(x,  z  =  rj,t)  is  the  vertical  velocity  at  the  surface.  Equations 
(23)  and  (24)  can  be  solved  if  an  associated  solution  of  Laplace’s  equation 
is  used  to  obtain  the  vertical  velocity.  Examples  of  successful  direct 
approaches  using  this  problem  formulation  may  be  found  in  [19],  [20]. 
These  approaches  often  utilize  methods  such  as  Fourier  transforms  to 
obtain  solutions  to  Laplace’s  equation  efficiently.  They  are  thus  usually 
restricted  to  periodic  domains  and  are  difficult  to  apply  to  the  general 
coastal  wave  propagation  problem. 

Alternately,  Radder  [21]  has  suggested  an  approach  in  which  the 
Hamiltonian  evolution  equations  are  written 

in  terms  of  rj  and  the  stream  function  at  the  free  surface, which  is  eval¬ 
uated  after  conformally  mapping  the  fluid  domain  into  a  strip.  Due  to 
the  transformation  method  utilized,  the  method  is  not  immediately  ap¬ 
plicable  to  wave  propagation  in  two  horizontal  directions.  Otta  &  Dinge- 
mans  [22]  show  a  method  for  computing  solutions  of  Radder’s  model 
utilizing  sine-series  to  discretize  the  remaining  Fourier  integral.  They 
further  use  a  check  on  the  approach  to  singularity  of  the  Jacobian  of  the 
transformation  from  physical  to  mapped  space  as  a  test  for  breaking, 
and  modify  the  surface  and  stream  function  evolution  subsequent  to  the 
break  point  to  model  the  effect  of  wave  breaking,  with  reasonable  success. 


Use  of  Hamilton’s  principle  to  construct  evolution  equations  in  the 
long  wave  setting  usually  proceeds  along  slightly  different  lines,  and  we 
defer  a  discussion  of  this  aspect  until  Section  4. 

2.3.2  Boundary  integral  equations. 

An  alternate  approach  to  the  inviscid  irrotational  flow  problem  results 
from  the  fact  that  Laplace’s  equation  for  the  fluid  interior  may  be  con¬ 
verted  to  an  integral  equation  on  the  surface  bounding  the  fluid  domain. 
Denoting  the  boundary  by  T,  which  is  either  a  closed  contour  in  the  case 
of  flow  in  one  horizontal  coordinate  and  z,  or  a  closed  surface  in  the  fully 
three-dimensional  case,  we  may  write  an  equation  for  the  value  of  the 
potential  at  a  point  on  the  boundary  surface  as 

<f>(xi)  =  jT  [<j>n(x)G(x,Xi)  -  <f>(x)Gn(x,  *,)]  dT(x)  (25) 

where  G  is  the  free-space  Green  function  and  where  subscript  n  denotes  a 
derivative  in  the  direction  of  the  outward  normal  vector.  Several  different 
approaches  exist  for  solving  (25)  together  with  the  surface  boundary  con¬ 
ditions.  Grilli  and  coworkers  [23]  [24]  use  the  dynamic  surface  boundary 
condition  (12)  together  with  a  kinematic  condition 

Dx  =  V3<£;  D()  =  ()t  +  V3  •  V3()  (26) 

to  step  the  solution  forward  in  time  and  update  the  position  of  the  bound¬ 
ary,  after  which  (25)  is  used  to  update  the  velocity  field.  Applications 
of  this  model  have  been  limited  to  date  to  two-dimensional  problems  in 
the  vertical  plane.  Grilli  et  al  [25]  have  compared  numerical  calcula¬ 
tions  of  solitary  waves  shoaling  on  a  plane  beach  to  laboratory  data,  and 
have  found  that  the  model  predicts  the  actual  wave  motion  to  within  the 
accuracy  of  the  available  data. 

Several  researchers  have  implemented  similar  approaches  using  this 
modelling  technique  [26]  [27]  [28]  [29].  Fully  three-dimensional  applica¬ 
tions  of  the  technique  are  typically  based  on  panel  methods  [30]  [31]  and 
are  still  relatively  scarce. 


2.4  Scaling  of  dimensional  variables 

Development  of  approximate  formulations  for  long  wave  motion  in  a  layer 
of  water  usually  rests  on  the  identification  of  two  scaling  parameters;  a 


water  depth  to  wavelength  ratio  (denoted  here  by  //),  and  a  wave  height  to 
water  depth  ratio  (denoted  here  by  e).  An  alternate  nonlinear  parameter 
can  be  obtained  by  multiplying  fi  and  e  to  obtain  the  ratio  of  wave  height 
to  wave  length,  the  so-called  wave  steepness. 

Following  [32] ,  we  use  a  reference  wavenumber  ko  to  scale  horizontal 
distances  x,  y,  a  reference  water  depth  h0  to  scale  the  vertical  coordinate  z 
and  local  depth  h(x,  y ),  and  amplitude  a  to  scale  the  surface  displacement 
tj.  We  then  introduce  the  parameters  8  —  a/ ho  and  y2  =  (koho)2.  Based 
on  these,  we  choose  a  scale  of  (Ar0 (flr^o)1/,2)— 1  for  time  t  and  8ho(gho)1^2  /  y 
for  velocity  potential  <j>.  Introducing  these  scales  into  the  boundary  value 
problem  for  inviscid,  irrotational  motion  leads  to  the  problem 


<t>zz  +  /i2VV  —  0; 

—h  <  z  <  8r] 

(27) 

<f>z  +  y2Vh  •  V<j)  =  0; 

z  =  —h 

(28) 

v  +  <t>t  +  ^[(W)2  +  ^(^)2] = 

<<) 

II 

(29) 

r, It  +  8V(f>  •  Vrj  -  <f>z  -  0; 

y2 

II 

^3 

(30) 

We  develop  an  equation  expressing  volume  flux  conservation  by  integrat¬ 
ing  (27)  over  z  from  —h  to  8rj  and  using  (28)  and  (30)  to  obtain 


rjt  +  V  •  M  =  0; 


(31) 


In  the  following,  we  use  (31)  to  obtain  expressions  for  mass  conservation, 
while  a  momentum  equation  is  obtained  using  the  Bernoulli  equation 
(29).  Choice  of  representative  length  and  time  scales  for  the  propagation 
of  irregular  waves  in  a  variable-depth  region  is  problematic,  since  waves 
satisfying  some  scaling  argument  or  restriction  in  one  portion  of  the  do¬ 
main  are  almost  certain  to  violate  these  conditions  in  another  portion 
of  the  domain.  We  proceed  here  with  a  standard  construction  of  the  di¬ 
mensionless  formulation.  We  will  address  the  effect  of  violating  the  basic 
scaling  assumptions  as  a  separate  issue  later. 


2.5  An  aside  on  nomenclature 

Approximate  models  for  wave  propagation  are  obtained  through  a  num¬ 
ber  of  routes,  and  confusion  may  arise  if  there  is  not  a  precise  under¬ 
standing  of  how  the  name  of  a  result  is  tied  to  the  method  for  obtaining 
that  result.  The  rules  followed  below  are  thus  outlined  here. 


Most  of  the  models  discussed  below  are  related  to  a  description  of  the 
flow  field  based  on  a  series  solution  of  Laplace’s  equation  retaining  terms 
to  0(p 2),  which  gives  a  vertical  distribution  of  horizontal  velocity  (at 
least  over  a  flat  bottom)  that  is  quadratic  in  2  or  some  related  vertical 
coordinate.  The  corresponding  vertical  distribution  of  vertical  velocity  is 
linear  in  z.  A  quick  check  of  the  approximate  expressions  for  the  velocity 
potential  used  below  shows  that  the  expression  for  the  fluid  vorticity 
is  identically  zero  to  0(p2),  and  thus  the  flow  fields  are  irrotational  to 
the  order  of  the  approximation.  We  will  refer  to  any  model  based  on 
a  potential  as  a  starting  point,  or  on  an  assumption  of  the  form  of  the 
internal  kinematics  that  would  correspond  to  the  use  of  a  potential,  as  a 
Boussinesq-type,  or  in  short,  a  Boussinesq  model.  Since  the  specification 
of  canonical  variables  in  the  Hamiltonian  formulation  uses  a  velocity 
potential,  approximate  models  obtained  using  this  approach  fall  in  this 
category  as  well. 

In  contrast,  there  are  a  number  of  models  that  assume  a  priori  that 
the  fluid  moves  in  vertical  columns  which  remain  vertical  and  are  at  most 
stretched  by  up  and  down  motion  of  the  surface.  This  assumption  has 
been  used  in  a  number  of  models,  including  models  based  on  approxima¬ 
tion  of  the  governing  Euler  equations  (Serre  [33],  Su  &  Gardner  [34])  or 
satisfaction  of  global  mass,  momentum  and  energy  constraints  (Green  & 
Naghdi  [35]).  These  models  are  variously  referred  to  as  Serre  equations 
or  GN  equations,  and  we  will  show  below  that  they  are  in  fact  the  same, 
at  least  at  the  leading  order  of  approximation  used  in  constructing  them. 

Finally,  we  will  make  a  distinction  between  fully-nonlinear  and  weakly- 
nonlinear  models.  We  will  refer  to  a  model  as  being  fully-nonlinear  if  the 
known  information  about  the  internal  flow  field  is  used  to  satisfy  the 
surface  boundary  conditions  imposed  at  the  free  surface,  with  the  only 
truncation  allowed  being  used  to  maintain  consistancy  in  the  ordering 
with  respect  to  the  parameter  p2.  In  this  context,  the  usual  Airy  or 
Nonlinear  Shallow  Water  (NSW)  equations  would  be  the  appropriate 
fully-nonlinear  equations  in  the  limit  p2  — >  0.  Any  model  which  invokes 
a  further  restriction  on  nonlinearity  by  fixing  a  relation  between  the  pa¬ 
rameters  8  and  p2  in  a  regime  where  p2  < Cl  will  be  referred  to  as  being 
weakly-nonlinear. 
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3  Derivation  of  Boussinesq  equations 


The  basic  problem  of  water  wave  propagation  is  two-dimensional  in  the 
horizontal  coordinates  x  =  (x,y),  and  we  will  start  with  models  that 
provide  a  full  treatment  of  arbitrary  motions  in  this  coordinate  system. 

The  central  goals  in  the  derivation  of  an  approximate  wave  propaga¬ 
tion  model  in  free  surface  hydrodynamics  are  to: 

1.  obtain  an  approximation  for  the  dependence  of  the  solution  on  the 
vertical  coordinate  z , 

2.  eliminate  the  cross-space  2  in  favor  of  a  model  in  the  propagation 
space  x,y,t. 

In  the  case  of  Boussinesq  models,  the  representation  of  the  vertical  de¬ 
pendence  of  the  solution  as  a  low-order  polynomial  leads  to  a  recursion 
relation  in  which  higher  order  coefficients  in  the  series  are  determined 
by  the  spatial  variation  of  the  lowest  order  coefficient.  The  choice  of 
this  coefficient  then  specifies  the  dependent  variable.  This  choice  is  non¬ 
unique,  and  thus  there  will  be  a  family  of  models  satisfying  a  given  set 
of  scaling  restrictions,  rather  than  one  fixed,  unique  Boussinesq  model. 
Mathematicians  have  typically  played  down  the  importance  of  the  dis¬ 
tinctions  between  the  members  of  this  family,  tending  to  lump  them 
together  as  representations  of  a  given  asymptotic  regime.  However,  as 
will  be  seen  below,  the  various  forms  of  the  model  equations  can  have 
significantly  different  properties  (particularly  in  terms  of  linear  disper¬ 
sion  characteristics)  when  applied  in  water  depths  which  are  large  rela¬ 
tive  to  the  asymptotic  range  of  validity.  For  this  reason,  modellers  are 
quite  interested  in  the  differences  in  behavior  of  the  various  members  of 
the  asymptotically-equivalent  family.  Witting  [36]  pointed  out  that  it 
is  possible  to  construct  weakly- dispersive  models  that  are  characterized 
by  dispersion  relations  having  the  form  of  rational  polynomial  or  Pade 
approximants,  and  that  are  numerically  far  more  accurate  as  the  depth 
becomes  relatively  large.  We  illustrate  the  consistant  construction  of  such 
a  model  below,  following  an  approach  developed  by  Nwogu  [37].  Other 
means  for  obtaining  models  with  extended  accuracy  in  their  dispersion 
relations  include  direct  manipulation  of  highest-order  (dispersive)  terms 
in  an  existing  set  of  governing  equations  [38]  [39],  or  the  use  of  an  al¬ 
ternate  vertical  shape  function  in  order  to  introduce  correct  dispersion 
properties  at  some  reference  frequency  [40]  [41]. 


The  derivation  of  the  Boussinesq  equations  illustrated  here  proceeds 
by  obtaining  an  approximate  solution  to  the  Laplace  equation  in  the  fluid 
interior  and  then  using  the  resulting  information  in  the  dynamic  surface 
boundary  condition  (29)  and  the  depth  integrated  continuity  equation 
(31).  As  pointed  out  by  Mei  [32],  the  resulting  derivation  does  not  require 
any  a  priori  assumption  about  the  size  of  all  information  provided  by 
(29)  and  (31)  can  be  retained.  The  resulting  fully-nonlinear  equations 
will  depend  on  the  level  of  accuracy  retained  in  the  approximation  of  (j). 
Weakly-nonlinear  Boussinesq  equations  are  then  obtained  by  invoking 
the  restriction  8  —  0(p2). 

3.1  Approximate  expressions  for  the  velocity  po¬ 
tential 

The  Boussinesq-type  equations  may  be  obtained  by  introducing  a  series 
expansion  for  (j)  of  the  form 

OO 

<Hx,z,i)  =  +  z)n<f>n(x,t)  (32) 

71=0 

An  expression  for  <f>  which  retains  terms  to  0(fi2)  and  satisfies  the  bottom 
boundary  condition  is  given  by 

<t>  =  4(x,0  -  fi2(h  +  z)Vh  •  V^o  -  ^(^-LfLvVo  +  0(p4)  (33) 

where  <j>o  is  the  value  of  the  velocity  potential  at  z  —  —h.  Various  forms 
of  the  evolution  equations  are  obtained  by  replacing  (f>0  by  the  value  of 
the  potential  at  any  level  in  the  water  column,  or  by  the  integral  over 
the  depth.  Any  choice  will  lead  to  a  set  of  model  equations  with  the 
same  level  of  asymptotic  approximation  but  with  numerically  different 
dispersion  properties,  as  discussed  below.  As  a  first  example,  we  will 
use  a  depth-averaged  value  of  the  potential  and  a  related  depth-averaged 
velocity,  following  Wu  [42]  and  Mei  [32].  We  define  the  depth-averaged 
potential  according  to 

^(x,  t)  =  jj  f_h  <Kx,  t)dz  (34) 

where  H  =  (h  +  Srj )  is  the  total  local  water  depth.  Using  (33)  in  (34) 
gives 


This  expression  is  then  used  in  (33)  to  obtain  an  expression  for  (f>  in  terms 
of 

2  2 

<f>==^>+lL.^H-2{h  +  z))Vh-V4>  +  ^(H2-3{h+zf)V24>  +  0{fJ,4)  (36) 

Alternately,  following  Wei  et  al  [43]  and  Chen  &  Liu  [44],  we  denote 
<j)a  as  the  value  of  <j>  at  a  reference  elevation  z  =  za(x),  or 

<j>a  =  <t>o~  »2{h  +  za)Vh  ■  V<t>0  -  +  0(n 4)  (37) 

This  expression  is  then  used  in  (33)  to  obtain  an  expression  for  <j>  in  terms 
of  (f>a . 

<t>  =<!>„  +  »2(z„  ~  *)V  •  (/>W„)  +  i/.2(4  -  z2)W0  +  0(fii)  (38) 


3.2  Model  based  on  depth- integrated  potential  or 
velocity 

3.2.1  Two-equation  model  for  rj  and  <j> 


Wu  [42]  was  one  of  the  first  to  point  out  that  a  two-equation  model  for  77 
and  the  chosen  dependent  variable  representing  the  potential  appeared  to 
be  a  good  candidate  for  computational  work,  since  it  reduces  the  number 
of  equations  and  dependent  variables  by  one  compared  to  a  formulation 
involving  r]  and  a  horizontal  velocity  vector.  Wu  formulated  a  model  in 
the  standard  Boussinesq  approximation,  including  the  effects  of  applied 
surface  pressure  and  moving  bottom.  We  illustrate  the  derivation  of  his 
model  here  and  provide  a  fully-nonlinear  version  as  an  intermediate  step. 


Using  (35)  in  (31)  to  obtain  an  expression  for  volume  flux  M  gives 


M  =  H 
+ 


W  +  V?  {\(^H  “  2VA)Vfc  •  V<£ 
^H{2VH  -  3V/i)V2^}]  +  0(n4) 


(39) 


The  expression  goes  to  zero  identically  as  the  total  depth  H  goes  to  zero, 
which  serves  as  a  natural  shoreline  boundary  condition.  The  correspond¬ 
ing  form  of  the  Bernoulli  equation  (29)  is 


lrD(Vh  ■  + 

z 


f  D(V2« 


Sfi2 

~6~ 


tf2(V2^)2  -  3 H(Vh  ■  V4>)V24>  -  3 (Vh  •  V^)2]  =  0{fi4)  (40) 


where  the  operator  D()  =  ()t  +  8V<j>-V().  The  usual  Boussinesq  approx¬ 
imation  is  recovered  by  retaining  terms  of  0(8 ,  fi2)  which  gives  the  set  of 
equations 

M  =  HVj>  -  fi2Vh  |  •  (hV$)  -  y  V2^|  +  0(82,  fy2, /)  (41) 


and 


•  (hVj,,)  - 


=  0(6\6li\y.')  (42) 


Equations  (41)  and  (42)  are  identical  to  the  results  given  by  Wu  [42]  after 
neglecting  pressure  forcing  and  bottom  motion. 


3.2.2  Three-equation  model  for  ??  and  u 

A  model  based  on  the  depth-averaged  horizontal  velocity  is  obtained 
next.  The  velocity  u  is  defined  by 

1  fSr> 

n  =  hLv*  (43) 

2  2 

=  V<£  +  7-(V#  -  2 Vh)(Vh  ■  V0)  +  ^(2 VH  -  ZVh)V24>  +  0(/x4) 
2  6 

Using  (43)  in  the  expression  (39)  for  the  mass  flux  gives 

M  =  Hu  +  0(n4)  (44) 

which  in  fact  is  accurate  to  all  orders  in  fi.  The  corresponding  momentum 
equation  is  obtained  by  taking  the  gradient  of  the  Bernoulli  equation  (40), 
which  gives 

u*  +  6( u  •  V)u  +  Vrj  -  [(VH  -  2 Vh) Vh  •  u]t 

|d-(V/.u)  +  ^D-(V-u) 

+^-V  [if2(V  •  u)2  -  3H(Vh  •  u)(V  •  u)  -  3(VA  •  u)2 

6  L 

-H(2VH  -  3 Vh)  ■  u(V  •  u)  -  3 (VH  -  2 Vh)  ■  u(Vh  •  u)]  =  0(/x4)( 45) 


[H(2VH  -  SVh)V  •  u]t  -  /i2V 


where  D*()  =  (),  +  8u  •  V()  =  DQ  +  0{fi2). 

The  weakly  nonlinear  Boussinesq  approximation  is  obtained  by  ne¬ 
glecting  terms  of  0(8  fi2)  or  higher  in  (45),  giving  (44)  for  the  mass  flux 
and 

u,  +  «(u- V)u  +  V,-/12  |v(V-(Mr,))-  jV(V-u<)  =0(6ti2,^1) 

(46) 

for  the  momentum  equations.  Equation  (46)  was  given  originally  by 
Peregrine  [9].  Retaining  full  nonlinear  effects  but  restricting  to  constant 
still  water  depth  h  =  h0  reduces  (45)  to 

2 

ut  +  8(u  •  V)u  +  Vt?  +  (tf2V2 u)t 

+/i2V  -^2u-V2u  +  ^2(V-u)2-^(V-u)t  =  0(n4)  (47) 

given  by  Mei  [32]. 


3.3  Model  based  on  potential  or  velocity  at  a  fixed 
elevation 


3.3.1  Two-equation  model  for  rj  and  <j>a 


Now,  we  derive  a  two-equation  model  for  rj  and  <j>a  following  Wei  et  al 
[43].  Using  (38),  the  expression  for 

M  in  (31)  becomes 


V2V<^  (48) 


M  =  H  V<^  +  /z2{V  zaV  ■  (hV<t>a)  +  -f' V2<f>, 


(h  -  Sri) 


V(V-(/>V^)) 


(, h 2  -  h&rj  ±  (Sri)2) 


The  corresponding  form  of  the  Bernoulli  equation  (29)  becomes 

V  +  <Aat  +  ^(V^)2  +  ^2[(^-^)V*(W^  +  ^(4-W2)VVat 
+  8 (I2  {V<f>a  •  [V^V  •  ( hV<f>a)  +  {za  -  8V)V(V  ■  (hV<f>a))] 

+  •  \zaVzaV2K  +  l-{z2a  -  (^)2)V(V2^) 

+  \  [V  •  (hV<t>a)}2  +  8VV  ■  (hV<j>a)V2<f>a  +  |(^)2(V2^)2}  =  0  (49) 
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Weakly-nonlinear  Boussinesq  equations  are  obtained  by  neglecting 
terms  of  0(8/x2).  The  modified  expression  for  volume  flux  M  is 


M  =  HV<f>a  +  fi2{hV 


•  (hv<pa)  +  -fV‘4>a 


+ yV(v  -ShVM)  -  j' ?2v  J 


(50) 


and  the  Bernoulli  equation  reduces  to 

r)  +  <f>at  +  f  [zav  •  (hv<i>at)  +  \zy2<f>at 


=  0 


(51) 


Equations  (50)  and  (51)  were  given  in  [44].  These  results  are  equivalent 
to  the  two-equation  model  (41)  -  (42)  of  Wu  [42]  to  within  rearrangements 
of  dispersive  terms. 


3.3.2  Three  equation  model  based  on  7/  and  ua 


We  further  introduce  a  horizontal  velocity  ua  as  ua  =  V<f> \Za.  Retaining 
terms  to  O(fP)  and  to  all  orders  in  8  gives  a  fully  nonlinear  version  of 
the  model  with  volume  flux 


M 


H 


Ua  +  1?  | -  ^( h 2  -  h8rj  +  (8rj)2)  V(V  •  ua) 

+  [*«  +  \{h  ~  8v)]  V(V  •  (feua))}]  +  0(m4)  (52) 

and  momentum  equation 

Uat  +  <5(u<*  •  V)u„  +  Wr]  +  /^2v !  +  8fi2\ 2  =  0(fl4)  (53) 


where 

Vx  =  i4V(V-u„,)  +  z«V(V -(/Hi.,)) 

-V  [i(5,)!V  •  u„,  +  ■  (Au„)] 

V2  =  V  [(*,  -  &,)(u„  •  V)(V  •  (feu„))  +  1(4  -  (6v)°)(ua  •  V)(V  ■  u.)' 
1 


The  weakly-nonlinear  Boussinesq  equations  of  Nwogu  are  recovered 

•  by  neglecting  terms  of  0(ix4,S(j,2),  yielding  the  expressions 

M  =  ifu0  +  ^{(^- jj  V(V  •  u.)  +  (hz„  +  j'j  V(V  •  (Au„))j 

(55) 

0  and 

uat+S(ua-V)ua+V7}+M2  |y  V(V  •  uQt)  +  zaV[V  •  (fcu*,)]}  =  0(<5//2^4) 

(56) 

•  The  fully  nonlinear  models  derived  here  all  have  mass  flux  M  — ►  0  at 
the  shoreline,  where  H  — *  0.  This  result  is  expected  on  physical  grounds 
and  appears  in  the  nonlinear  shallow  water  equations  and  in  Boussinesq 
models  where  the  depth-averaged  velocity  is  the  dependent  variable.  This 
condition  is  not  automatically  satisfied  by  Nwogu’s  or  any  weakly  nonlin- 

•  ear  Boussinesq  model  based  on  a  velocity  other  than  the  depth-averaged 
value,  making  the  application  of  these  models  problematic  at  the  shore¬ 
line.  All  fully-nonlinear  variations  of  any  of  the  possible  model  systems 
should  recover  this  condition  correctly. 


3.4  Choice  of  za  and  linearized  dispersion  relations 

As  the  limitation  to  relatively  shallow  water  /it2  < C  1  is  the  primary 
factor  affecting  applicability  of  the  Boussinesq  models,  it  is  important 
to  understand  the  impact  that  the  choice  of  dependent  variables  has  on 
the  implied  linear  dispersion  relation  associated  with  each  model.  To 
investigate  this,  we  consider  linearized  forms  of  the  model  based  on  the 
depth- averaged  potential  (41)- (42)  and  the  potential  chosen  at  a  fixed 
depth  za ,  (50)-(51).  For  the  case  of  one- dimensional  propagation  in  water 
of  constant  depth  h  =  ho,  the  former  may  be  written  as 


7/t  “1“  $xx  -  6 

(57) 

"b  1)  ixt  ~  0 

(58) 

and  has  the  associated  dispersion  relation 

2  1 

C  1  +  I/*’ 

(59) 
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where  c  is  the  phase  speed  normalized  by  the  nondispersive  phase  speed 
\fgh.  Chen  h  Liu’s  model  (50),  (51)  is  written  in  1-D,  linearized  form  as 


Vt  +  4*  a  xx  +  H2(a  +  -)4>  axxxx 

(frat  T  f]  "b  ft  &4>axxt 

and  has  the  corresponding  dispersion  relation 


0 

0 


c2  =  1  -  (Q  +  Dm2 


1  —  a//2 


(60) 

(61) 


(62) 


In  (60)-(62),  a  is  given  by 


1  (za\2  .  za 
°  2\h)  +  h 


(63) 


The  appearance  of  the  free  parameter  a  in  (62)  gives  the  model 
described  here  (or  any  other  extended- accuracy  Boussinesq  model)  the 
ability  to  predict  accurate  phase  speeds  over  a  wide  range  of  water  depths. 
For  example,  the  choice  a  =  —1/3  (corresponding  to  za  =  (l/\/3  —  1  )ft  = 
—0.423ft)  reproduces  the  relation  (59)  based  on  the  depth-averaged  ve¬ 
locity.  As  shown  originally  by  Witting  [36],  a  better  result  is  obtained 
by  forcing  the  approximate  formula  to  coincide  with  the  (2,2)  Pade  ap- 
proximant,  given  by  the  choice  a  =  —2/5  (which,  in  turn,  corresponds 
to  za  =  (l/\/5  —  1  )h  =  —0.553 h).  Finally,  Nwogu  [37]  pointed  out  that 
more  accurate  results  can  be  obtained  at  intermediate  depth  by  choos¬ 
ing  a  such  that  some  measure  of  error  is  minimized.  Nwogu  chose  to 
minimize  the  mean  square  error  between  (62)  and  the  exact  result 


c2  = 


tanh  n 


over  the  range  0  <  u>2h/g  <  k,  and  obtained  the  result 


a  =  —0.39; 


-0.53ft 


(64) 


(65) 


The  dispersion  relations  resulting  from  use  of  the  depth-averaged  velocity, 
the  Pade  approximant  and  Nwogu’s  choice  of  o;  are  compared  to  the  exact 
result  in  Figure  1.  Further  refinements  have  been  suggested  by  Chen  & 
Liu  [44],  who  minimized  the  mean  square  error  in  a  combination  of  the 
phase  and  group  velocities  and  obtained  a  slightly  different  value  of  a. 
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Figure  1:  Dependence  of  linear  phase  speed  c  on  dimensionless  relative 
depth  //.  Solid  line,  linear  theory  (64).  Dashed  line,  Boussinesq  equa¬ 
tions  based  on  depth-averaged  velocity  (59).  Dash-dot  line,  Boussinesq 
equations  with  a  =  —2/5  in  (62).  Dotted  line,  Boussinesq  equations  with 
a  =  -0.39  in  (62). 


4  Boussinesq  approximations:  Hamiltonian 
formulations 


As  an  alternative  to  the  standard  derivation  outlined  above,  several  au¬ 
thors  have  provided  derivations  of  weakly-nonlinear  Boussinesq  equa¬ 
tions  based  on  the  Hamiltonian  and  resulting  canonical  evolution  equa¬ 
tions.  These  results  are  reviewed  here  in  the  simplified  context  of  two- 
dimensional  propagation  in  water  of  constant  depth  h.  A  detailed  review 
of  existing  variable  depth  formulations  may  be  found  in  Dingemans  [45]. 

Most  of  the  work  on  Boussinesq  approximations  through  Hamiltonian 
formulations  follows  along  lines  laid  out  by  Broer  and  associates  [46]  [47]  [48]. 
The  Hamiltonian  may  be  written  as  the  sum  of  kinetic  and  potential  en- 
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ergy  contributions 


H  =  T  +  V  (66) 

Introducing  the  non-dimensionalization  of  section  2.4  and  normalizing  Ti 
by  pga2  / 2,  we  get  (in  dimensionless  form) 


T 

V 


=  //xKv^+>>! 

=  fin'd* 


dzdx 


(67) 

(68) 


The  difficulty  in  deriving  the  canonical  evolution  equations  lies  in  speci¬ 
fying  the  canonical  momentum  £(x,  t)  in  terms  of  the  velocity  potential 
4>{x,  z  =  r),  t).  Broer  and  most  subsequent  researchers  have  proceeded  by 
splitting  the  kinetic  energy  integral  into  an  integral  from  the  bottom  to 
the  mean  water  level  plus  an  integral  from  the  mean  water  level  to  r/: 


T  =  %  +  T„ 


(69) 


where 


%  = 


Tr, 


dzdx 


/iXKv*+>>: 

=  //X  Kv*+>> 

=  «/ /^[vfv^fc+oiV/) 


dzdx 


(70) 

(71) 


by  virtue  of  the  fact  that  <j>z  =  0(g2)  at  leading  order,  and  where  we 
introduce  a  Taylor  expansion  about  z  =  0  and  retain  terms  appropriate 
to  a  weakly-nonlinear  formulation.  (Note  that  this  approach  appears  to 
eliminate  the  possibility  of  obtaining  a  fully-nonlinear  formulation,  as 
given  in  section  3).  Subsequent  development  depends  on  the  assumption 
that  the  value  of  the  potential  at  the  surface  may  be  represented  in  terms 
of  the  value  at  the  still  water  level  z  =  0.  Letting  £o(x,  t )  =  <^(x,  z  =  0,  t), 
we  may  use  Green’s  identity  together  with  Laplace’s  equation  and  the 
bottom  boundary  condition  to  replace  (70)  by 


%  =  ^J  Jj0(h),=odx  (72) 

It  is  illustrative  to  consider  the  implications  of  the  solution  to  the  full 
linear  problem  in  the  development  of  the  weakly-dispersive  approxima¬ 
tion  here.  For  waves  propagating  in  constant  depth  in  two  dimensions, 
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the  full  solution  to  the  problem  in  an  unbounded  domain  may  be  written 
(in  dimensional  form) 


«x’ f  >  = ; hJL  <73> 

where  k  =  |k|  is  the  magnitude  of  the  wavenumber  vector  and  £0  is  the 
Fourier  transform  of  the  potential  at  the  still  water  level.  The  linearized 
kinetic  energy  may  then  be  written  as  [45]  [50] 

To  =  |£o|2ktanh  khdk  (74) 

Alternately,  using  Parseval’s  theorem,  we  can  write  (74)  in  terms  of  an 
integral  over  space,  as 

% = \  I L  ve°  ■  dx  ^ 

where  R/  is  a  real,  positive-definite  operator  with  a  Fourier  transform  (or 
symbol)  given  by 

-  tanh  kh 


For  the  case  of  weak  dispersion  (tt  <  1)  we  can  introduce  the  approxi¬ 
mation 

tanh  kh  ■  1,,,,,  ,__n 

— =  1  -  3 (kh)2  +  O(kh)4  +  ...  (77) 

which  corresponds  to  an  operator  R  of  the  form 

R  =  1  +  \h2V2  (78) 

which  is  not  a  positive-definite  approximation  to  the  full  integro-differential 
operator  by  virtue  of  the  fact  that  (77)  takes  on  negative  values  at  large 
kh.  The  central  problem  in  developing  the  Boussinesq  approximation  is 
then  to  obtain  a  replacement  for  R  of  suitable  form,  so  as  to  maintain 
positive-definiteness  as  well  as  to  obtain  relatively  accurate  dispersive 
properties. 

We  now  return  to  nondimensional  variables  and  evaluate  (70)  directly 
and  compare  to  the  results  obtained  above.  Choosing  za  =  0  as  the 
reference  level  in  (38),  we  retain  terms  to  0(fi 4)  in  (f>  and  obtain  the 


results  (neglecting  bottom  slope) 
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(f>(x.,z,t)  =  &  -  fi2(hz  +  y )V2£o 

V  (|jA4  -  j h\h  +  zf  +  ±(h  +  zy)  V2V%  (79) 


Also  note  that  £  =  <f>(z  =  r})  =  £0  +  0(6n2),  so  that  £0  may  be  replaced 
by  £  in  the  Hamiltonian  for  the  lowest  order  Boussinesq  model.  Using 
(79)  in  (72)  then  gives 


T0  =  J  Jvt-ihRVOdx  (80) 

where 

R=l+/i2jV2;  R=l-^2  (81) 

in  agreement  with  the  first  two  terms  in  the  expansion  (77).  Combining 
(68),  (71)  and  (80)  gives  the  expression 


•  (6r)  +  fcR)V£dx 


(82) 


Use  of  (21)  and  (22)  then  yields  the  canonical  evolution  equations 

<h  =  -V-{HV(]-SjV2V2(,  (83) 

(,  =  jVf.Vf-1,  (84) 

The  problems  imposed  by  the  non-positive  definiteness  of  the  operator 

(81)  and  resulting  Hamiltonian  (82)  may  be  seen  by  examining  the  dis¬ 
persion  relation  for  the  linearized  version  of  (83)-(84),  given  by 

c2  =  1  -  (85) 

Note  that  c  takes  on  imaginary  values  for  values  of  //  =  kh  >  y/Z,  indi¬ 
cating  the  instability  of  short  wave  solutions  to  the  problem.  This  is  a 
large  value  of  p  by  the  standards  originally  applied  to  Boussinesq-type 
models  but  falls  well  within  the  range  of  accurate  prediction  capabil¬ 
ity  of  models  with  extended  accuracy.  In  addition,  since  short  waves 
on  the  scale  of  2Ax  are  certain  to  arise  in  numerical  calculations,  the 
present  model  may  not  be  used  as  the  basis  for  numerical  calculations 


unless  heavy  short-wave  damping  is  artificially  imposed.  This  damping, 
in  turn,  makes  the  replication  of  steep  nonlinear  wave  features  impossible 
in  a  practical  sense. 

It  is  thus  necessary  to  alter  (82)  in  such  a  way  that  short  wave  in¬ 
stabilities  are  avoided.  From  an  ad-hoc  point  of  view,  this  end  may  be 
achieved  simply  by  replacing  the  operator  R  by  another  having  a  more 
suitable  symbol,  and  hence  a  well  behaved  linear  dispersion  relation;  for 
example,  the  choice 


Ri  = 


_  1  -  (a  +  1/3  V  _ 


Ri  = 


1  +  (a  +  l/3)h3V2 
1  +  ah3V2 


(86) 


1  —  afj, 2 

reproduces  the  extended  dispersion  relation  of  Nwogu’s  model  (62)  and 
gives  the  modified  evolution  equations 

(1  +  n2aV2)Vt  =  -V  •  [HV£]  -  fx2(a  +  i)fc3V2V2£  +  0(8  fi2,^)  (87) 

together  with  (84).  Most  authors  take  the  additional  step  of  replacing 
(82)  with  a  quadratic  form,  thus  guaranteeing  positive-definiteness.  For 
example,  Mooiman  [50]  introduces  the  quadratic  form 


H  =  I  j  [H(GV02  +  l2}  dx 

which  corresponds  to  the  canonical  evolution  equations 

m  =  -v  •  (G'(ffG(V«) 


(88) 

(89) 

(90) 


where  G*  is  defined  such  that  G*G  =  R/  to  the  required  order  of  ap¬ 
proximation.  Mooiman  makes  the  choice 


G  —  Gq>6  —  Lfc  1L0; 


1  -  /A2V2 
p  6 


(91) 


Requiring  G  to  reproduce  the  leading  order  Taylor  series  behavior  of  Ri 
in  (81)  leads  to  the  condition 


b  —  a  —  1 

and  a  resulting  linear  dispersion  relation 

I ±mv 


c‘  =  G‘  = 


l  +  l/*2. 


(92) 


(93) 
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Based  on  comparisons  with  the  exact  dispersion  relation,  Mooiman  sug¬ 
gested  the  choice  of  parameters  b  =  1.9,  a  =  0.9.  However,  as  shown  in 
Figure  2,  it  is  not  clear  that  this  is  the  best  choice  over  a  range  of  values 
of  n  =  kh  appropriate  to  the  intermediate  depth  problem. 

Values  closer  to  6  =  1.7,  a  =  0.7  would  seem  to  provide  a  better  fit.  It 
is  also  clear  that  models  that  reproduce  the  Pade  form  of  the  dispersion 
relation  provide  a  more  accurate  representation  of  the  dispersion  relation 
over  this  same  range. 


Figure  2:  Squared  linear  phase  speed  for  model  equations  relative  to 
squared  phase  speed  for  exact  relation  (64).  Solid  line,  Pade  approximant 
(62)  with  a  =  —2/5.  Dashed  line,  Mooiman  model  (93)  with  b  =  1.9. 
Dash-dot  line,  Mooiman  model  (93)  with  b  =  1.8.  Dotted  line,  Mooiman 
model  (93)  with  b  =  1.7. 

Additional  derivations  of  weakly-nonlinear  Boussinesq  models  have 
been  provided  by  Van  der  Veen  &  Wubs  [51]  and  Yoon  &  Liu  [52].  In 
particular,  the  model  of  Van  der  Veen  k,  Wubs  reproduces  the  Pade  form 
of  the  linear  dispersion  relation  and  may  thus  be  expected  to  predict  re- 


suits  with  comparable  accuracy  to  those  of  the  Nwogu  or  other  extended 
weakly-nonlinear  model  equations. 


5  Serre  equations 

An  independent  line  of  approach  to  the  derivation  of  evolution  equations 
has  been  pursued  by  several  authors  starting  with  Serre  [33].  In  this 
approach,  a  model  for  the  internal  flow  field  is  adopted  in  which  the 
horizontal  velocity  is  assumed  to  be  uniform  over  depth  and  the  vertical 
velocity  varies  linearly  from  the  bottom  to  the  free  surface.  Models  are 
then  developed  by  obtaining  appropriate  integrals  over  depth  of  the  Euler 
equations  and  applying  surface  boundary  conditions.  The  dependent 
velocity  variable  is  thus  the  depth-averaged  velocity  u,  which  coincides 
with  the  local  horizontal  velocity  over  the  entire  water  column.  To  date, 
applications  of  the  resulting  approach  have  been  limited  to  one  horizontal 
dimension  (taken  to  be  x).  Extensive  details  on  model  derivations  may 
be  found  in  Dingemans  [45]. 

Note  that  the  assumed  form  of  the  kinematics  in  these  models  implies 
that  the  computed  flow  field  is  rotational.  For  motion  in  the  x,  z  plane, 
we  have  vorticity  u>  =  dw/dx  =  0(p2),  and  thus  vorticity  enters  the  flow 
field  description  at  the  same  order  as  the  leading  order  dispersive  terms. 

Serre-type  models  are  derived  assuming  a  particular  relation  between 
scaling  parameters  8  and  fi.  For  the  case  where  the  usual  weakly-nonlinear 
Boussinesq  approximation  is  employed  (8  =  0(fi2)),  Dingemans  [45] 
shows  that  the  derivation  reproduces  the  equations  of  Peregrine  [9]  writ¬ 
ten  in  terms  of  depth-averaged  velocity.  Serre  [33]  and  Su  &  Gardner  [34] 
employed  the  relation  8  =  O(^),  which  leads  to  the  retention  of  some  non¬ 
linear  effects  in  dispersive  terms  in  the  resulting  model  equations.  For 
the  case  of  a  horizontal  bottom,  Serre  and  Su  &  Gardner  obtained  the 
set  of  equations  [54] 

Tjt  +  ( Hu)x  =  0  (94) 

2  1 

ut  +  8uux  +  tjx  +  -fi2HxT  +  -fi2HTx  =  0  (95) 

o  o 

where  u  is  the  x-component  of  u  and 

T  =  H(8(u2x  -  uuxx)  -  uxt )  (96) 
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Equations  (94)- (96)  have  an  exact  solitary  wave  solution,  given  by 


with  K  =  (3<$/4(l  +  8))1!2  and  c  —  1/(1  +  8 ) 1^2  and  where  the  nonlin¬ 
ear  parameter  8  describes  the  ratio  of  wave  height  over  water  depth.  In 
comparison,  K  for  the  solitary  wave  solution  of  the  KdV  equation  (132) 
is  given  by  K  —  3<$/4.  The  Serre  solution  produces  solitary  wave  crests 
which  are  quite  broad  relative  to  either  exact  solutions  or  to  solutions 
based  on  the  weakly  or  fully-nonlinear  Boussinesq  models  of  the  previous 
section.  In  Figure  3,  we  show  a  comparison  of  solitary  waves  based  on 
the  accurate  numerical  solution  of  Tanaka  [53]  (solid  lines),  the  weakly 
nonlinear  KdV  solution  (dashed  line),  a  numerical  solution  of  the  fully 
nonlinear  Boussinesq  equations  (52)  -  (53)  [55]  (dash-dot  line),  and  the 
present  solution  of  the  Serre  equation  (dotted  line).  The  figure  shows 
that  the  fully-nonlinear  Boussinesq  model  provides  a  highly  accurate  re¬ 
production  of  the  exact  solitary  wave  shape  up  to  very  high  wave  heights. 
In  contrast,  the  wave  crest  predicted  by  Serre  equations  is  quite  broad 
and  becomes  broader  as  height  increases,  relative  to  any  other  solution 
technique. 

Seabra-Santos  et  al  [54]  have  provided  a  variable-depth  form  of  the 
Serre  equations  given  by 


H(ut  +  8uux)  + 


rjt  +  (Hu)x 


0  (99) 

Hhx(  i  +  s  +  irxioo) 


with 

E  —  hx^zif  -j-  8uux)  8hxxu  vlOl) 


Finally,  we  note  that  the  linearized  dispersion  relation  correspond¬ 
ing  to  the  models  described  here  is  the  relation  (59),  which  requires 
improvement  before  the  models  can  be  used  to  model  waves  in  interme¬ 
diate  depths.  Steps  in  this  direction  have  been  made  [66]  but  are  largely 
undocumented  to  date. 


Figure  3:  Comparison  of  solitary  wave  shapes  for  wave  heights  8  = 
0.4, 0.6, 0.8.  Solid  line  -  Tanaka’s  solution.  Dashed  line  -  KdV  soli¬ 
tary  wave  solution.  Dash-dot  line  -  numerical  solution  of  fully  nonlinear 
Boussinesq  equations  (52)  -  (53).  Dotted  line  -  Serre  equation  solution 
(97)  -  (98) 

6  Green-Naghdi  (GN)  equations 

The  final  method  for  obtaining  approximate  governing  equations  consid¬ 
ered  here  is  the  method  of  Green  &  Naghdi  [56]  [35].  In  this  approach,  an 
assumption  is  made  about  the  kinematic  properties  of  the  velocity  field, 
after  which  conservation  laws  governing  the  coefficients  of  the  velocity 
field  are  derived.  The  original  approach  to  the  problem  [56]  used  the 
theory  of  directed  fluid  sheets  or  Cosserat  surfaces.  Green  &  Naghdi  [35] 
provide  a  parallel  derivation  of  their  approximate  equations  based  on  the 
explicit  conservation  laws  for  inviscid  fluid  flow.  These  equations  are 
now  commonly  referred  to  as  Green-Naghdi  or  GN  equations.  In  this 
formulation,  the  method  takes  on  something  of  the  flavor  of  a  Galerkin 
approximation  of  the  full  problem,  using  the  assumed  form  of  the  velocity 
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field  as  the  shape  function.  The  resulting  models  are  fully  nonlinear,  as 
no  scaling  assumptions  are  made  about  the  relative  height  to  depth  ratio 
and  all  boundary  data  are  evaluated  at  the  instantaneous  free  surface. 

The  original  GN  model  uses  the  same  kinematic  assumption  as  the 
Serre  model;  a  linear  variation  of  vertical  velocity  over  depth  and  a  depth- 
uniform  horizontal  velocity.  The  equations  of  Green  &  Naghdi  can  be 
written  in  the  present  notation  following  Miles  &  Salmon  [57];  see  also  [58] 

%  +  V-{Hu)  =  0  (102) 

u* +  <5u- Vu  + Vt/  =  —pi2 A  (103) 

where 

A  =  j^(H2D2v)  (104) 

for  uniform  depth  and 

A  =  i  {V  [HD2(2r)  -  ft)]  +  {Vrj)D2( 2rj  -  h)  +  ( Wh)D2{2h  -  V)} 

6  ^  '  (105) 

for  variable  depth,  where  D()  =  ()t  +  u  •  V()  and  where  we  retain  the 
notation  u  for  the  uniform-over-depth  horizontal  velocity.  As  is  the  case 
with  the  Serre  equations,  (102)-(103)  represents  a  system  having  non¬ 
zero  vorticity.  Miles  &  Salmon  [57]  consider  the  vorticity  and  show  that 
a  potential  vorticity  may  be  derived  and  is  conserved  by  the  system 
for  any  motion  starting  from  rest.  Calculations  for  upstream-advancing 
solitons  generated  by  ships  moving  at  transcritical  speed  have  been  made 
by  Ertekin  et  al  [58].  There  has  been  very  little  direct  comparison  to 
data  for  this  particular  model  system. 

Owing  to  the  similarity  in  the  assumed  velocity  profiles,  one  would 
suspect  that  the  Serre  and  GN  models  should  bear  a  striking  resemblance. 
To  investigate  this,  we  may  rewrite  (102)  as 


Dr)  =  —Hux  (106) 

for  the  case  of  constant  depth  and  propagation  in  one  dimension.  The 
expression  for  D2r]  then  becomes  (in  1-D) 

D2rj  =  H[8(u2x  —  uuxx)  —  uxt]  =  T  (107) 

where  T  is  defined  in  (96).  Expanding  the  one-dimensional  version  of 
(103)  then  gives  (95)  identically,  proving  the  equivalence  of  the  models. 
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Based  on  this  correspondence,  we  may  write  a  two-dimensional  version 
of  T  as 

r  =  H  {$[( V  •  u)2  -  u  •  V(V  •  u)]  -  V  •  ut}  (108) 

thus  extending  the  Serre  equations  to  two  horizontal  dimensions. 

As  is  the  case  with  the  other  weakly  dispersive  models  described  here, 
the  leading  order  GN  model  needs  to  be  extended  to  include  more  ac¬ 
curate  linear  dispersion  effects  before  it  can  be  applied  to  intermediate 
depth  propagation.  Shields  &  Webster  [59]  [60]  have  provided  a  frame¬ 
work  for  such  an  extension.  They  proceed  by  writing  the  velocity  field 
(in  the  present  notation)  as 


K 

u(x,  *,<)=£  W„(x,  t)sn  (109) 

n=0 

where  K  is  a  finite  integer  corresponding  to  the  level  of  the  resulting 
model,  and  where  s  is  the  vertical  position  in  a  mapped  coordinate  sys¬ 
tem  which  places  the  bottom  at  s  =  —  1  and  the  water  surface  at  s  =  1. 
Equations  are  then  found  which  govern  the  evolution  of  the  Wn,  and 
which  satisfy  mass  conservation  and  the  kinematic  surface  boundary  con¬ 
dition.  The  Euler  equations  are  satisfied  in  an  approximate  sense,  using 
a  weak  variational  formulation  in  which  the  basis  functions  of  (109)  are 
used  as  weighting  functions.  The  resulting  system  consists  of  K  equa¬ 
tions  for  the  amplitudes  of  the  velocity  field  together  with  an  equation 
governing  the  change  in  thickness  of  the  water  column.  Further  details 
may  also  be  found  in  Demirbilek  &  Webster  [61].  To  date,  these  theories 
have  only  been  implemented  for  one  horizontal  direction. 

Shields  &  Webster  [59]  have  shown  calculations  of  solitary  waves  and 
regular  cnoidal  waves  for  GN  models  up  to  level  three,  and  have  shown 
that  convergence  towards  numerically  exact  results  for  wave  shape  and 
wave  speed  is  more  rapid  than  for  corresponding  perturbation  series. 
Calculations  of  waves  shoaling  in  variable  depth  are  also  provided  by 
Shields  [62],  who  shows  good  agreement  between  wave  height  evolution 
and  wave  shape  for  shoaling  periodic  waves,  in  comparison  to  laboratory 
data  [63].  Additional  example  calculations  are  provided  by  Demirbilek 
&  Webster  [64].  Webster  &  Wehausen  [65]  have  also  applied  the  method 
to  the  problem  of  resonant  Bragg  reflection  of  surface  waves  by  undular 
bed  features. 


7  Model  performance 


The  testing  of  existing  modelling  schemes  has  been  fragmentary,  espe¬ 
cially  with  regard  to  model  intercomparisons  relative  to  single  data  sets 
or  more  exact  calculations.  The  principal  exception  to  this  is  provided  by 
Dingemans  [66],  who  has  provided  a  comparison  of  calculations  of  nine 
models  to  a  single  laboratory  data  set.  The  main  issues  to  be  addressed 
in  testing  of  the  models  include  accuracy  of  intermediate  depth  linear 
propagation  (including  dispersion,  shoaling,  refraction  and  diffraction), 
and  reproduction  of  nonlinear  wave  features  during  the  final  stages  of 
shoaling  and  wave  breaking.  In  addition,  the  question  of  whether  ex¬ 
tended  Boussinesq  models  give  accurate  predictions  of  nonlinear  features 
in  intermediate  water  depths  is  of  some  importance,  and  has  not  been  ad¬ 
equately  addressed  to  date.  An  outline  of  available  results  and  significant 
findings  is  provided  here. 


7.1  Intermediate  depth  propagation 

Several  experiments  have  been  performed  which  provide  accurate  data 
sets  for  intermediate  depth  waves  propagating  over  bathymetry.  The  data 
set  of  Berkhoff  et  al  [67],  for  wave  focussing  over  a  submerged  shoal, 
has  been  extensively  used  to  test  the  performance  of  parabolic  model 
schemes  [68].  A  schematic  of  the  bottom  geometry  is  shown  in  Figure 
4.  Kirby  and  Dalrymple  [69]  showed  that  the  experimental  results  are 
significantly  affected  by  nonlinearity,  and  that  parabolic  model  schemes 
that  account  for  the  leading  order,  self-interaction  effect  in  Stokes  wave 
theory  are  capable  of  reproducing  the  focusing  behavior  (including  the 
effect  of  nonlinear  defocussing)  observed  in  the  experiment.  In  addition, 
waves  in  the  experiment  propagate  from  a  relatively  deep  region,  where 
kh  ~  2,  to  much  shallower  depths.  This  experiment  is  thus  ideal  for 
testing  the  intermediate  depth  shoaling  and  nonlinear  properties  of  the 
extended  Boussinesq  schemes  developed  above.  Wei  and  Kirby  [70]  per¬ 
formed  such  a  comparison  using  the  weakly-nonlinear  Nwogu  equations. 
Their  published  results  were  inaccurate,  due  to  problems  with  the  inci¬ 
dent  wave  boundary  condition  as  well  as  a  misinterpretation  of  whether 
the  model  had  reached  an  asymptotic  steady  periodicity.  Corrected  re¬ 
sults  for  a  comparison  of  measured  and  computed  wave  heights  for  the 
measurement  transects  shown  in  Figure  4  are  shown  in  Figure  5.  The 


results  show  close  agreement  with  data,  with  prediction  accuracy  similar 
to  that  available  using  large- angle  parabolic  model  schemes  [68].  These 
results  are  not  changed  in  any  significant  way  when  the  fully  nonlinear 
formulation  is  used.  Similar  results  using  a  model  based  on  Hamilton’s 
principal  have  been  shown  by  Mooiman  [71]. 


Figure  4:  Geometry  of  submerged  shoal  experiment  of  Berkhoff  et  al  [67] 
including  measurement  transects  (from  [70]). 


7.2  Solitary  waves  on  a  plane  beach 

Waves  close  to  the  break  point  can  often  have  height  to  depth  ratios  larger 
than  one  during  their  rapid  unsteady  evolution,  and  thus  the  restrictions 
of  a  weakly  nonlinear  theory  are  dramatically  violated  during  the  final 
stages  of  shoaling.  The  problem  of  wave  evolution  close  to  breaking  thus 
provides  an  excellent  test  of  the  realtive  accuracy  of  weakly-nonlinear  and 
fully-nonlinear  versions  of  the  Boussinesq  models.  Wei  et  al  [43]  have 
studied  the  shoaling  of  solitary  waves  up  to  the  break  point,  and  have 
compared  results  from  the  weakly-nonlinear  model  (55)  -  (56)  and  the 
fully-nonlinear  model  (52)  -  (53)  to  results  obtained  using  an  accurate 
boundary  integral  solution  [24].  The  solitary  wave  test  is  a  good  choice 


Figure  5:  Comparison  of  measured  (circles)  and  computed  crest-to- 
trough  wave  heights  for  submerged  shoal  experiment  [67]  [70]. 

for  performing  model  comparisons,  since  the  wave  can  be  started  inside 
the  model  domain  and  there  are  no  initial  transients  or  other  features 
in  the  solution  resulting  from  boundary  treatments,  which  can  cloud  an 
evaluation  of  the  basic  models. 

Wei  et  al  considered  solitary  waves  of  three  initial  heights,  each  shoal¬ 
ing  over  4  different  beach  slopes.  Results  showed  that  weakly-nonlinear 
Boussinesq  calculation  consistantly  showed  an  overprediction  of  wave 
height  as  waves  approached  the  break  point,  together  with  overly  narrow 
and  curved  crests  (and  resulting  high  horizontal  velocities.)  In  contrast, 
the  fully-nonlinear  model  equations  provided  an  accurate  prediction  of 
wave  height  evolution  (see  Figure  6)  as  well  as  accurate  reproduction  of 
wave  shape  and  magnitude  of  maximum  horizontal  velocities.  The  pre¬ 
dicted  point  where  horizontal  velocity  at  the  wave  crest  coincides  with 
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Figure  6:  Comparison  between  boundary  integral  model  (solid  lines), 
weakly  nonlinear  Boussinesq  model  (dashed  lines)  and  fully  nonlinear 
Boussinesq  model  (dash-dot  lines)  of  shoaling  rates  H/h  for  solitary 
waves  with  8  =  (A)  0.20,  (B)  0.40,  (C)  0.60  in  (a),(b),(d)  and  8  =  (A) 
0.30,  (B)  0.45,  (C)  0.60  in  (c);  shoaling  on  a  slope:  (a)  1:100;  (b)  1:35;  (c) 
1:15;  (d)  1:8.  Circle  symbols  denote  locations  of  the  breaking  point  for 
which  the  wave  has  a  vertical  tangent  on  the  front  face  in  the  boundary 
integral  computations.  (From  [43]). 
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the  crest  translational  speed,  which  would  be  an  indication  of  limiting 
wave  height  from  a  physical  point  of  view,  agreed  closely  with  the  pre¬ 
dicted  point  of  the  onset  of  breaking  in  the  boundary  integral  calcula¬ 
tions.  Overall,  the  significance  of  the  additional  nonlinear  effects  in  the 
final  stages  of  shoaling  was  quite  clear. 

7.3  Additional  model  intercomparisons 

Dingemans  [66]  has  provided  a  comparison  between  the  predictions  of 
seven  weakly  dispersive  wave  models,  a  boundary  integral  model,  a  fully- 
dispersive  model  based  on  canonical  Hamiltonian  evolution  equations, 
and  a  laboratory  data  set.  The  lab  experiment  considered  the  shoaling 
of  an  initially  sinusoidal  regular  wave  over  a  bar  with  sloping  sides.  The 
bar  height  occupied  75%  of  the  water  depth  at  the  crest.  Waves  shoaled 
significantly  over  the  bar  and  broke  over  the  bar  crest  in  one  of  the  three 
cases  studied.  The  strongly  deformed  waves  then  propagated  back  into 
relatively  deep  water  on  the  lee  side  of  the  bar,  and  the  waves  decomposed 
into  a  fundamental  and  free  higher  harmonics  with  only  weak  nonlinear 
interactions.  The  experimental  data  thus  provides  a  strong  test  of  the 
linear  dispersive  characteristics  of  the  models  as  well  as  their  ability  to 
reproduce  strong  nonlinear  distortions  during  the  shoaling  process.  Tests 
with  breaking  were  not  a  primary  focus  of  the  comparison  study.  For  a 
case  with  moderately  long  initial  waves  (2s  period  in  0.8m  of  water),  the 
study  showed  fairly  conclusively  that  the  primary  consideration  affecting 
model  performance  was  the  accuracy  of  the  model  dispersive  properties, 
with  the  retention  of  higher  order  nonlinear  effects  being  a  distinguishing 
feature  only  after  correction  of  the  dispersive  properties.  In  this  instance, 
the  best  model  performance  in  the  group  of  long  wave  models  was  found 
in  an  undocumented  Serre  model  with  corrected  dispersive  properties  due 
to  Barthelemy  &  Guiborg  of  LEGI-IMG  (see  [66]  for  details  of  model 
formulation).  This  conclusion  unfortunately  did  not  survive  to  the  next 
case  studied,  which  used  a  shorter  but  steeper  incident  wave.  In  this 
instance,  the  weakly-nonlinear  Boussinesq  models  of  Madsen  et  al  [39] 
and  Mooiman  [50]  outperformed  the  Serre  model  with  corrected  linear 
dispersion.  This  result  was  unexpected  but  may  be  another  indication 
of  the  Serre-type  model  tending  to  deviate  from  accurate  solutions  as 
steepness  increases,  as  indicated  in  the  solitary  wave  comparison  above. 
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8  Additional  physical  effects 


The  models  described  to  this  point  represent  the  propagation  of  waves 
seaward  of  the  surf  zone,  in  the  range  where  friction  is  negligible  and 
energy  is  essentially  conserved.  A  comprehensive  model  of  nearshore  pro¬ 
cesses  needs  to  be  able  to  describe  the  additional  complexity  of  surf  zone 
wave  behavior,  including  wave  breaking  and  runup  on  the  shoreface.  The 
breaking  process  is  not  described  by  the  basic  framework  given  above. 
In  addition,  the  process  of  fluid  motion  in  the  swash  zone  occurs  in  a 
regime  where  the  dispersive  properties  of  the  Boussinesq  model  are  not 
appropriate.  Nevertheless,  there  has  been  a  great  deal  of  progress  made 
in  extending  the  range  of  weakly  dispersive  models  into  the  surfzone, 
and  thus  providing  a  comprehensive  model  of  wave  motion  in  the  coastal 
environment. 

8.1  Wave  breaking 

The  application  of  the  Boussinesq  models  in  the  surf  zone  region  is  prob¬ 
lematic  from  the  viewpoint  of  consistency.  Nonlinearity  becomes  strong 
(6  — >  0(1))  while  dispersive  effects  should  disappear  (yu2  — *■  0).  In  this 
limit,  the  Boussinesq  model  would  approach  the  usual  NSW  equations. 
These  equations  predict  the  eventual  steepening  and  “breaking”  (mani¬ 
fested  through  the  convergence  of  solution  characteristics)  of  any  initial 
wave  form,  and  thus  do  not  allow  for  the  modelling  of  non-breaking  waves 
over  any  appreciable  distance.  The  NSW  equations  provide  a  success¬ 
ful  and  well  tested  framework  for  modelling  the  dissipation  and  runup 
of  surfzone  waves  [72]  [73]  [74]  [75].  The  numerical  approach  is  usu¬ 
ally  based  on  the  Lax-Wendroff  scheme  or  similar  dissipative  schemes 
which  preserve  the  mass  and  momentum  conserving  properties  of  the 
governing  equations,  but  which  dissipate  energy.  This  approach  pro¬ 
vides  predictions  of  surf  zone  wave  heights,  fluid  velocities  and  skewness 
and  asymmetry  statistics  which  are  in  reasonable  agreement  with  mea¬ 
surements  [76].  The  approach  has  the  disadvantage  that  the  modeller 
does  not  have  control  over  the  dissipative  processes  involved,  which,  in 
turn,  are  not  specified  by  any  physically  based  model.  Finally,  the  NSW 
equations  cannot  accurately  predict  the  propagation  of  incident  waves 
seaward  of  the  surfzone. 

In  contrast,  the  schemes  utilized  to  solve  the  Boussinesq  equations 


35 


typically  are  at  least  intended  to  be  conservative  of  mass,  momentum 
and  energy,  and  do  not  form  weak  solutions  or  shocks  as  the  wave  fronts 
steepen.  Instead,  the  increased  wave  front  curvature  leads  to  an  increase 
in  the  effect  of  dispersive  terms  as  water  depth  decreases,  rather  than 
a  decrease.  Since  the  schemes  are  basically  conservative,  the  equations 
respond  to  the  increased  curvature  by  radiating  shorter  waves  which  are 
slower  and  hence  trail  behind  the  main  wave  crest.  This  effect  is  equiv¬ 
alent  to  the  undular  bore  phenomenon  seen  in  low  Froude  number  hy¬ 
draulic  jumps  [77],  the  difference  being  that  breaking  would  not  occur 
and  replace  the  tendency  to  radiate  the  short  waves  as  wave  height  in¬ 
creased  beyond  some  value.  It  is  up  to  the  modeller  to  invoke  a  dissipation 
mechanism  and  include  it  in  the  model  equations. 

A  number  of  different  models  for  wave  breaking  effects  have  been  used 
to  date.  These  models  can  be  roughly  separated  into  two  groups:  models 
where  the  rate  of  dissipation  depends  on  some  geometric  or  global  (in 
time  or  space)  property  of  the  wave  train,  and  models  where  the  rate 
of  dissipation  is  determined  entirely  by  local  (in  time  or  space)  proper¬ 
ties.  Within  each  group,  it  is  further  possible  to  distinguish  two  basic 
techniques  for  achieving  a  damping  effect:  use  of  an  eddy  viscosity  for¬ 
mulation,  or  use  of  an  applied  pressure  distribution. 

Eddy  viscosity  models  have  the  longest  history  in  application.  These 
involve  extending  the  momentum  equation  by  the  addition  of  a  dissipa¬ 
tion  term,  indicated  schematically  by 

ut  +  8uux  +  T)x  —  ( vbux)x  +  dispersive  terms  =  0  (HO) 

Note  that  some  authors  [78]  [79]  write  the  additional  term  with  the  eddy 
viscosity  coefficient  outside  the  second  derivative.  This  implies  that  a 
contribution  to  flow  momentum  is  imposed  when  breaking  occurs,  in 
contrast  to  the  momentum-conserving  bore  process  in  the  non-dispersive 
theory.  In  situations  where  dissipation  is  imposed  globally  and  spatial 
variations  in  vb  over  a  wave  period  are  small,  this  effect  is  minor.  How¬ 
ever,  at  the  onset  of  breaking  or  in  models  where  dissipation  is  localized 
and  spatial  gradients  of  viscosity  are  large,  this  momentum  source  effect 
can  be  quite  severe,  and  should  be  avoided  by  correctly  specifying  the 
dissipation  term. 

Karambas  et  al  [78]  provide  an  example  formulation  where  vb  is 
calculated  based  on  a  mixing  length  hypothesis  and  then  applied  globally 
over  the  waveform.  No  comparisons  to  data  are  provided. 
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Zelt  [80]  has  utilized  a  spatially  localized  eddy  viscosity  formula¬ 
tion  [81]  based  on  a  mixing  length  hypothesis.  Following  Zelt,  vb  is  given 
by 

vb  =  -l2ux  (111) 

with  the  mixing  length  l  given  by 

l  =  B'rH  (112) 

Heitner  and  Housner  [81]  chose  7  =  2  based  on  comparisons  between 
numerical  results  and  experiments  on  the  width  of  the  bore  region  in  a 
steady  hydraulic  jump.  The  factor  B  is  used  to  introduce  a  breaking 
criterion  and  is  given  by 
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where  u*  is  a  critical  velocity  gradient  taken  to  be 
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In  practice,  the  expression  (112)  is  linearized  to  give  /  =  B^h,  and  a 
smoothing  filter  is  used  to  smooth  the  computed  vb  distribution  slightly. 
Zelt  [80]  used  this  formulation  in  the  context  of  a  weakly- nonlinear  Boussi- 
nesq  model  in  Lagrangian  coordinates  and  studied  the  breaking  and 
runup  of  solitary  waves  on  a  plane  slope,  and  obtained  reasonable  agree¬ 
ment  with  laboratory  data  [82].  More  recently,  Wei  and  Kirby  [55]  have 
used  this  method  to  study  the  breaking  of  a  random  wave  train  over 
a  plane  slope.  Wei  and  Kirby  used  the  laboratory  data  of  Mase  and 
Kirby  [83],  who  considered  the  propagation  of  a  wave  train  corresponding 
to  a  Pierson-Moskowitz  spectrum  over  a  1:20  beach  slope.  An  example 
of  30  seconds  of  run  time  is  shown  in  Figure  7,  where  the  initial  mea¬ 
surement  at  a  depth  of  47cm  is  shown  together  with  shoaled,  breaking 
waves  at  h  —  17.5,15,12.5,10,7.5,5  and  2.5cm  progressing  up  the  fig¬ 
ure.  These  tests  were  conducted  with  a  peak  frequency  of  about  1  Hz, 
corresponding  to  a  value  of  kh  ~  2  in  the  deeper  part  of  the  tank.  The 
Nwogu-type  model  with  improved  dispersion  was  used  as  the  basis  for 
calculations.  Computations  using  the  model  based  on  depth-averaged 
velocity  (44)  and  (46)  were  also  attempted  and  failed  to  reproduce  the 
shoaling  and  propagation  of  the  shorter  wave  components  in  the  wave 
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train,  as  would  be  expected  for  waves  starting  in  this  large  a  depth.  The 
calculations  shown  here  reproduce  the  arrival  time  of  individual  waves  as 
well  as  the  height  and  general  shape  of  waves  in  the  surfzone. 

Statistical  measures  based  on  the  entire  experimental  series  (covering 
about  800  waves)  were  computed.  Figure  8  shows  computed  and  mea¬ 
sured  third-moment  statistics  for  an  entire  run,  and  indicates  that  the 
simple  eddy  viscosity  formulation  is  capable  of  providing  accurate  pre¬ 
dictions  of  wave  skewness  and  asymmetry,  which  are  important  for  the 
subsequent  calculation  of  sediment  transport  properties. 


o 
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Figure  7:  Comparison  of  surface  elevation  at  different  water  depths.  Solid 
line  -  experimental  data  [83].  Dashed  line  -  numerical  calculation  based 
on  (52)  -  (53).  Isolated  blips  on  records  under  wave  records  indicate 
temporally  localized  breaking.  (From  [55]). 

Karambas  &  Koutitas  [79]  have  provided  a  much  more  elaborate 
scheme  for  computing  the  eddy  viscosity  v\,  using  an  energetic  eddy  length 
scaling  rather  than  a  mixing  length  based  on  the  total  depth.  They  obtain 
good  prediction  of  decaying  wave  heights  in  the  surfzone,  as  compared 


Figure  8:  Comparison  of  skewness  (circles)  and  -asymmetry  (stars)  at 
different  water  depths.  Solid  lines  -  experimental  data  [83].  Dashed  lines 
-  numerical  computations  based  on  (52)  -  (53).  (From  [55]). 


to  regular  wave  data  from  several  sources.  They  also  provide  setup  cal¬ 
culations  which  show  a  large  underprediction  relative  to  measured  data. 
This  underprediction  may  well  be  due  to  the  incorrect  positioning  of  the 
eddy  viscosity  in  the  model  dissipation  term  in  their  formulation,  which 
would  lead  to  a  modification  of  the  cross-shore  time-averaged  momentum 
balance. 


An  alternate  formulation  based  on  a  surface  roller  concept  [84]  [85] 
has  been  described  by  Schaffer  et  al  [86].  Using  a  weakly- nonlinear 
Boussinesq  formulation  based  on  surface  displacement  and  depth-integrated 
volume  flux,  they  formulate  the  model  (written  in  the  present  notation) 
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where  P  =  Hu  and  where  R  is  a  pressure- like  imposed  force  specified  by 
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the  roller  model.  R  is  found  to  be  an  excess  momentum  effect  which  can 
be  written  as  2 

R  =  M (117) 

where  M  is  the  actual  momentum  flux  computed  using  the  true  verti¬ 
cal  distribution  of  horizontal  velocity,  including  breaking  and  turbulence 
effects.  The  main  assumption  in  the  model  development  thus  revolves 
around  the  choice  of  the  assumed  velocity  profile.  Schaffer  et  al  choose 
a  simple  two-layer  model  where  the  roller  region  is  taken  to  be  a  simple 
quiescent  mass  of  water  translating  at  the  phase  speed  of  the  underlying 
wave.  Note  that  in  this  formulation,  where  R  takes  on  non-zero  values 
only  near  the  front  of  breaking  wave  crests,  the  overall  average  momen¬ 
tum  balance  is  not  affected  by  the  additional  term. 

Numerical  calculations  with  this  model  also  provide  good  predictions 
of  wave  height  decay,  and  additionally  provide  accurate  prediction  of 
surfzone  setup.  Additional  results  using  this  formulation  may  be  found 
in  Madsen  et  al  [87].  The  same  general  approach  has  also  been  used  by 
Brocchini  et  al  [88]  in  both  a  Boussinesq  and  Serre  model. 

8.2  Runup 

Modelling  the  complete  wave  transformation  process  in  the  cross-shore 
direction  requires  a  treatment  of  the  movement  of  the  water  line  on  the 
inclined  shore  face.  The  problem  of  modelling  runup  on  the  shore  face 
does  not  introduce  major  new  physical  considerations  in  the  modelling 
scheme  (unless  percolation  effects  are  to  be  considered).  Instead,  the 
problems  presented  are  primarily  numerical  and  are  associated  with  the 
movement  of  the  boundary  of  the  fluid  domain  in  a  fixed  horizontal  co¬ 
ordinate  system.  There  are  several  approaches  to  the  solution  of  this 
problem  which  are  described  in  the  literature,  which  can  be  categorized 
roughly  as 

1.  coordinate  transformation  techniques,  in  which  a  mapped  coordi¬ 
nate  follows  the  water  edge 

2.  grid  draining  and  filling  techniques 

3.  techniques  which  treat  the  entire  modelled  domain  as  part  of  the 
active  fluid  domain. 


The  first  category  of  models  is  exemplified  by  models  developed  in 
Lagrangian  coordinates.  For  the  particular  case  of  purely  cross-shore  mo¬ 
tion,  use  of  Lagrangian  coordinates  renders  the  shoreline  stationary,  and 
thus  the  runup  problem  is  essentially  solved.  Pedersen  and  Gjevik  [89] 
reported  calculations  of  the  runup  of  solitary  waves  using  a  Boussinesq 
model  written  in  Lagrangian  coordinates,  and  obtained  good  agreement 
with  laboratory  data.  This  approach  was  extended  by  Zelt  [80],  who 
incorporated  the  breaking  model  described  above  and  studied  the  runup 
of  both  breaking  and  non-breaking  solitary  waves,  again  obtaining  good 
results.  This  approach  can,  in  principal,  be  extended  to  two-dimensional 
problems  having  longshore  as  well  as  cross- shore  shoreline  movement. 
Use  of  a  purely  Lagrangian  scheme  for  surfzone  applications  is  problem¬ 
atic,  however,  since  the  mean  flow  developed  by  breaking  waves  would 
tend  to  advect  the  grid  longshore  relative  to  the  stationary  region  of 
interest.  It  is  possible  to  adopt  semi-Lagrangian  schemes,  where  only 
the  cross-shore  motion  is  treated  in  a  Lagrangian  framework.  Alter¬ 
natively,  simple  coordinate  stretching  techniques  may  be  utilized  which 
strain  the  model  grid  in  the  cross-shore  direction  in  order  to  follow  the 
moving  shoreline.  This  approach  has  been  used  by  Ozkan  &  Kirby  [90] 
in  schemes  for  longshore  runup  using  the  NSW  equations. 

The  second  category  of  shoreline  treatments  has  a  long  history  in  the 
context  of  storm  surge  and  tsunami  runup  calculations.  In  this  approach, 
the  water  level  in  the  most  shoreward  active  grid  in  the  numerical  fluid 
domain  is  checked.  The  shoreline  is  moved  landward  if  the  water  level 
exceeds  a  critical  value  or  the  grid  is  drained  and  inactivated  if  the  water 
level  falls  below  a  second  critical  value.  This  technique  is  most  applicable 
in  explicit  numerical  schemes  which  are  not  impacted  by  the  change  of  a 
grid  dimension.  Liu  et  al  [92]  provide  a  recent  example  of  this  approach, 
again  in  the  context  of  long  wave  runup  in  the  NSW  equations. 

Finally,  there  have  been  several  approaches  aimed  at  maintaining 
an  active  computational  grid  over  an  entire  pre-specified  model  domain 
containing  subaerial  regions.  Madsen  et  al  [87]  describe  the  application 
of  a  method  due  to  Tao  [93]  in  which  the  solid  seabed  is  replaced  by 
a  network  of  very  narrow  vertical  slots  superimposed  on  the  original 
topography.  Water  is  assumed  to  occupy  all  of  the  network  of  slots  up 
to  the  original  still  water  level,  and  thus  the  entire  model  grid  is  active 
from  the  start.  As  the  shoreline  advances,  slots  fill  until  the  water  level 
reaches  the  original  bed  level,  after  which  water  flows  over  the  newly 
submerged  bed.  In  a  sense,  the  method  is  analogous  to  coupling  the 


fluid  domain  to  a  porous  bed,  in  which  the  level  of  porosity  is  then  made 
as  small  as  possible.  Madsen  et  al  show  a  comparison  for  regular  wave 
runup  to  the  analytic  solution  of  Carrier  &  Greenspan  [91]  and  show 
that  the  level  of  apparent  porosity  in  the  numerical  scheme  must  be  very 
carefully  controlled  in  order  to  obtain  accurate  runup  calculations.  They 
also  found  that  it  was  neccessary  to  turn  off  the  dispersive  terms  in  the 
model  as  the  runup  tip  is  approached. 

An  alternate  approach  to  the  problem  of  maintaining  an  active  fluid 
domain  everywhere  uses  the  approach  of  imposing  a  bottom  friction  term 
which  becomes  extremely  large  as  the  local  total  depth  approaches  zero, 
thereby  “freezing”  a  thin  layer  of  fluid  on  an  otherwise  subaerial  re¬ 
gion  [94]. 

Following  this  approach,  Wei  &  Kirby  [55]  added  a  bottom  friction 
term  of  the  form 

F  -  u“lUal 

Fb  -  -~cpr  (118) 

and  maintained  an  active  grid  over  the  entire  domain  by  imposing  a 
minimum  total  depth.  This  method  also  requires  a  suppression  of  the 
model  dispersive  terms  at  the  runup  tip,  since  the  large  curvature  of  the 
water  surface  as  the  fluid  jumps  from  near-zero  to  finite  depth  tends  to 
generate  large  short-wave  oscillations  on  the  scale  of  the  grid  spacing. 


8.3  Waves  on  currents  and  mean  current  genera¬ 
tion 

Although  the  Boussinesq  model  is  most  often  thought  of  as  a  wave  prop¬ 
agation  model,  it  contains  the  framework  for  the  computation  of  steady 
or  low-frequency  motions  driven  either  by  imposed  currents  at  bound¬ 
aries  or  by  the  time-average  of  products  of  short  wave  components.  This 
may  be  seen  directly  by  splitting  the  velocities  and  surface  displacement 
into  wavy  and  slowly  varying  components,  and  then  averaging  the  result¬ 
ing  equations  to  obtain  separate  equations  for  waves  and  currents.  As 
an  illustration  of  the  results,  consider  the  decomposition  of  the  surface 
displacement  T]  and  velocity  field  u  into  a  current-induced  part  and  a 
wave-induced  part 


where  we  will  use  tensor  notation  for  convenience.  Peregrine’s  [9]  equa¬ 
tions  may  be  written  as 

rjt  +  [Hu,].  =  0  (120) 

ui,t  +  8xijUitj  +  ?7,i  —  n  =  0  (121) 

We  further  introduce  the  notation 

H  =  h  +  8r)  =  h  +  8{t]c  +  tjw )  —  Hc  +  8t]w  (122) 

Using  (119)  and  (122)  in  (120)-(121)  and  using  a  suitable  average  over 
wave  phase  denoted  by  (),  we  obtain  an  equation  for  mass  conservation 
for  the  mean  flow  given  by 


HCjt  +  (Hcucj  +  8(riwuwj))tj  —  0  (123) 


where  the  wave  averaged  quantity  is  the  wave-induced  mass  flux.  Denot¬ 
ing  the  total  mean  flow  velocity  by 


Utt  —  Ttc{  T  ( TJwUwi ) 

■tir 


we  may  further  write  the  mean  flow  momentum  equation  as 


(Hcuti)t  T  {HcutiUtj),j  -f-  Hcrjc<i  T  8S%jtj  —  0(8  ) 


(124) 


(125) 


where  dispersive  effects  on  the  slowly  varying  mean  flow  are  simply 
dropped  in  lieu  of  a  more  thorough  scaling  argument.  The  tensor  5,-j 
is  the  radiation  stress  tensor 

Sij  =  Hc(uwiuwj)  +  ^( r)l )  (126) 

and  the  term  Sijj  provides  the  forcing  to  drive  steady  or  low-frequency 
motions.  Turning  to  the  equations  for  the  fluctuating  motion,  we  re¬ 
turn  to  vector  notation  and  denote  the  wave-induced  velocity  as  u„  and 
current  velocity  as  uc.  The  resulting  governing  equations  are 

Vw,t  +  V -{Huw  +  8r]wuc- 8(r)wuw)]  =  0  (127) 

u**  +  <5(uc  +  u^)  •  Vuw  +  <$11^  •  Vuc  +  Vr/W 

-\^h2V(V-uwt)  =  0  (128) 

O 
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where  we  have  further  neglected  bottom  slope  effects  in  the  dispersive 
terms.  The  decomposition  here  is  based  on  the  assumption  that  the 
largest  component  of  the  mean  flow  is  of  the  same  order  as  the  wave 
induced  velocity;  =  0(8y/gh)  in  dimensional  form.  However,  it  is  oc¬ 
casionally  the  case  (particularly  in  tidal  inlets  and  estuarine  entrances) 
that  the  mean  flow  velocity  can  become  the  same  order  of  magnitude  as 
the  free  wave  speed  y/gh ,  and  thus  exceed  the  wave-induced  orbital  mo¬ 
tion  in  size.  This  case  is  not  adequately  covered  by  the  weakly-nonlinear 
Boussinesq  models  as  given  because  it  implies  a  fluid  velocity  which 
is  larger  than  assumed  in  the  initial  perturbation  analysis.  Yoon  and 
Liu  [95]  have  considered  this  case  in  detail  and  have  rederived  the  depth- 
averaged  equations  including  a  separate  large  current  from  the  start.  The 
principle  revisions  to  the  resulting  governing  equations  occur  only  in  the 
dispersive  terms  for  the  wave  motion,  and  essentially  occur  due  to  the 
advection  of  the  dispersive  wave  motion  by  the  large  current.  Equation 
(128)  may  be  corrected  to  this  more  general  level  of  approximation  by 
replacing  the  0(/x2)  dispersive  term  by  the  expression 


\g2H2cV  [V  •  uwt  +  uc  •  V(V  •  u.)]  (129) 

which  essentially  replaces  the  local  derivative  in  the  original  expression 
with  a  total  derivative  following  the  mean  flow.  The  fully  nonlinear 
Boussinesq  models  described  in  section  3  do  not  impose  an  absolute  scal¬ 
ing  on  the  size  of  velocity  components,  and  thus  contain  Yoon  &  Liu’s 
extension  of  the  theory  within  the  single  set  of  model  equations.  These 
models  are  thus  completely  valid  for  the  computation  of  wave  motions  on 
relatively  large  currents.  The  consistency  of  the  model  equation  with  the 
dispersive  term  (129)  may  be  seen  by  constructing  its  linear  dispersion 
relation,  which  reads  (in  dimensional  form) 


,  _  gVHc 
i  + 1  (kHcy 


(130) 


where 


cr  =  u  —  kuc 


(131) 


is  the  relative  or  intrinsic  frequency.  Equation  (130)  is  a  consistant  ap¬ 
proximation  to  the  full  linear  dispersion  relation  for  waves  on  a  large 
current  uc. 

Since  the  forcing  for  steady  or  low-frequency  motion  implied  by  (125) 
is  embedded  in  the  original  Boussinesq  or  Serre  model,  computations  in¬ 
volving  wave  breaking  should  generate  longshore  currents  or  other  more 


complex  nearshore  circulation  patterns  directly.  For  the  case  of  a  periodic 
incident  wave  train,  the  steady  wave-induced  current  can  be  extracted 
using  a  simple  time  average  of  the  computational  results.  An  example  of 
such  a  current  pattern  may  be  found  in  Sprensen  et  al  [96],  who  illus¬ 
trate  the  computation  of  a  rip  current  driven  by  longshore  bathymetric 
variations  on  an  otherwise  straight  coastline.  A  computed  water  surface 
and  one  half  of  the  resulting  circulation  cell  are  shown  in  Figures  9  and 
10. 


Figure  9:  Birds  view  of  the  free  surface  elevation  for  the  rip  channel  case 
with  regular  waves.  The  surface  rollers  are  shown  in  white,  (from  [96]). 

9  Models  for  wave  propagation  with  a  prin¬ 
cipal  direction 

In  addition  to  the  models  discussed  above,  which  are  all  essentially  fully 
two-dimensional  in  the  horizontal  plane  (or  have  the  capability  to  be 
extended  in  such  a  way),  a  considerable  amount  of  work  has  been  done 
on  modelling  schemes  which  involve  the  a  priori  choice  of  a  preferred 
or  principal  propagation  direction.  These  modelling  schemes  typically 
take  the  form  of  first-order  wave  equations  in  this  principal  direction,  ex¬ 
tended  to  include  effects  due  to  nonlinearity,  frequency  dispersion,  weak 


transverse  structure  in  the  wave  field  and  variations  in  the  fluid  domain 
in  the  propagation  or  transverse  directions.  This  subject  has  also  been 
recently  reviewed  by  Akylas  [97]. 

Historically,  the  first  member  of  this  class  of  models  was  developed 
by  Korteweg  and  deVries  [3]  and  is  known  as  the  Korteweg- deVries  or 
KdV  equation.  For  waves  propagating  in  the  + re-direction,  the  equation 
may  be  written  in  stationary  coordinates  as 

3  1 

Vt  +  Vx  +  ^7;Tlrlx  + I^2tVxxx  =  o  (132) 

2  o 

where  the  constant  depth  h  =  ho  and  where  the  linearized  phase  speed  c 
normalized  by  Co  =  y/gh 0  is  given  by 

e  =  1  -  i/i2  (133) 

Note  that  the  KdV  equation  exhibits  the  same  kind  of  cutoff  in  prediction 
of  a  positive  phase  speed  at  high  frequencies  as  seen  in  the  Boussinesq 
equations  based  on  velocity  at  the  mean  water  level,  although  the  prob¬ 
lem  of  imaginary  speeds  and  subsequent  instability  is  avoided.  This  effect 
was  addressed  by  Peregrine  [77]  and  Benjamin  et  al  [98],  who  introduced 
the  model  equation 

3  1 

Vt  +  Vx  +  =  o  (134) 

which  has  a  corresponding  linear  phase  speed 

_  1 

1  +  |  g2 

which  is  positive  for  all  wave  lengths.  Equation  (134)  is  commonly  re¬ 
ferred  to  as  the  BBM  or  Regularized  Long  Wave  (RLW)  equation. 

For  the  case  where  a  transverse  y-direction  structure  in  the  wavefield 
is  permitted  (with  characteristic  y  direction  scales  of  0(1/ fx)  relative 
to  re-direction  scales),  Kadomtsev  and  Petviashvili  [99]  introduced  the 
perturbed  KdV  equation  (now  commonly  referred  to  as  the  KP  equation) 

3  1  1 

(Vt  +  Vx  +  b-^Wx  +  g2^rixxx)x  +  l*2-jpYY  =  0  (136) 

where  Y  =  fxy.  As  in  the  case  of  the  KdV  equation,  the  KP  equation  is 
also  solvable  by  means  of  inverse  scattering  transforms  and  is  thus  the 


(135) 


Figure  10:  Depth-integrated  velocity  for  the  situation  shown  in  Figure  9 
focussing  on  a  circulation  cell  (from  [96]). 

subject  of  considerable  mathematical  interest.  There  is  also  a  solution 
describing  families  of  biperiodic  waves  (i.e.,  periodic  in  both  the  x  and  y 
directions);  see  Segur  and  Finkel  [100]  and  references  therein.  Hammack 
et  al  [101]  [102]  have  studied  the  correspondence  between  these  solutions 
and  short  crested  nonlinear  waves  generated  in  a  laboratory  basin,  and 
have  shown  a  close  correspondence  between  model  and  data,  and  have 
shown  that  the  theoretical  restriction  that  the  transverse  lengthscale  be 
large  relative  to  the  lengthscale  in  the  propagation  direction  is  not  a 
strong  restriction  on  applicability  of  the  model  equation. 

In  the  following,  we  review  recent  efforts  made  in  the  direction  of 
extending  and  applying  these  one- dimensional  or  weakly  two-dimensional 
models  to  the  problem  of  wave  propagation  in  regions  with  varying  depth. 
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9.1  Weakly  two-dimensional  propagation 

The  first  KP  equation  for  wave  propagation  in  a  variable  domain  with  ar¬ 
bitrary  depth  /i(x)  was  given  by  Liu  et  al  [103],  who  obtained  it  primarily 
from  an  examination  of  existing  KP  and  variable-coefficient  KdV  models. 
Using  a  more  formal  perturbation  expansion,  Philip  [104]  suggested  that 
a  more  appropriate  form  of  the  model  can  be  written  as 

(v t  +  cr)x  +  (c!V)y  =  0  (137) 

where  c  is  dimensionless  phase  speed  c  =  (h/ho)1^2,  and  where  we  as¬ 
sume  that  bottom  slope  in  the  propagation  direction  is  0(8).  This  equa¬ 
tion  is  the  same  as  the  equation  presented  in  [103]  to  within  terms  of 
0(fi2)  times  the  bottom  slope,  which  is  smaller  than  terms  retained  in 
the  model  equation.  The  model  was  used  as  the  basis  for  constructing  a 
parabolic  approximation  for  forward-scattered  waves  [68],  and  was  found 
to  be  an  accurate  predictor  of  wave  focussing  over  variable  bathymetry. 
Kirby  et  al  [105]  used  this  model  to  compute  the  Mach  reflection  of 
solitary  waves  from  a  vertical  wall,  and  further  extended  the  model  co- 
efficents  to  include  the  case  of  two-layer  flow  over  topography.  The  form 
of  the  linearized,  nondispersive  model  corresponding  to  (137)  can  also 
be  justified  using  a  splitting  method  applied  to  the  usual  linear  long 
wave  equation  [105]  [106];  this  approach  also  provides  the  equation  for 
the  backward-propagating  wave  component  as  well  as  the  coupling  be¬ 
tween  the  forward  and  backward  modes  due  to  linear  reflection  effects. 
A  similar  KP  model  has  been  developed  by  Chen  [107]. 

A  number  of  additional  studies  have  provided  KP  models  for  a  wide 
range  of  physical  settings,  including  surface  and  internal  waves  with 
Coriolis-induced  rotation  [108]  [109],  gradually- varying  channel  width 
and  background  flow  [110]  [111].  Chen  &  Liu  [113]  have  provided  a 
derivation  of  a  comprehensive  KP  model  covering  all  of  these  effects, 
with  some  restriction  on  the  allowed  variability  of  bottom  topography. 

In  most  of  the  studies  mentioned  above,  the  domain  in  which  waves 
are  propagating  is  assumed  to  be  a  channel  with  topography  which  varies 
very  weakly  in  the  transverse  direction.  The  depth  may  then  be  written 
as 

fe(x)  =  h(x )  +  8B(x)  (138) 

where  it  is  assumed  that  deviations  from  the  average  depth  in  the  trans¬ 
verse  direction  are  0(8)  relative  to  the  average  depth.  This  assumption 
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prompts  a  reformulation  of  the  bottom  boundary  condition  (6),  which 
after  expanding  about  the  depth  h  becomes 

uhx  +  •  VB  =  —  w  +  6Bwz  +  0(82)]  z  =  —h(x)  (139) 


This  expansion  has  the  effect  of  eliminating  transverse  depth  variation 
effects  in  the  leading  order  solutions  for  the  dependent  variables  in  the 
problem.  The  primary  consequence  of  this  is  to  remove  the  depth  h  from 
within  the  y  derivative  terms  in  (136)  or  (137).  The  resulting  model  for 
surface  waves  without  background  current  [112]  [113]  may  be  written  in 
the  present  notation  as 


'  „  o~hx  3  c6  2ch2  ccB  ^ 

Vt  +CT]X  +  fi  C~~~Tj  +  -~-T}T]x  +  V  -T~Vxxx  + 
t  4  h  la  o  lh  j 


where  c  =  (h/ho)1^2. 


For  simplicity,  most  studies  to  date  have  provided  computational  ex¬ 
amples  based  on  the  assumption  that  channel  sidewalls  are  impermeable 
vertical  structures.  Mathew  h  Akylas  [114]  provide  an  extension  to  the 
model  system  to  account  for  sidewalls  which  are  sloped  in  such  a  way  as 
to  still  be  confined  to  a  region  which  is  narrow  compared  to  the  main 
region  of  the  channel,  but  which  are  sloped  gradually  enough  to  affect 
the  dynamics  of  the  motion  in  the  channel. 


Published  tests  of  these  model  equations  against  actual  data  are  lack¬ 
ing.  The  model  (140)  with  the  assumption  of  weak  transverse  dependence 
in  the  bottom  topography  has  the  advantage  of  being  transformable  into 
a  standard  KP  equations  for  several  topographies  of  interest  [112]),  and 
thus  analytic  information  on  the  deformation  of  solitary  waves  by  this 
class  of  topographies  is  immediately  available.  It  is  the  present  author’s 
feeling  that  models  which  invoke  this  restriction  may  be  too  restrictive  for 
application  to  general  coastal  topographies,  although  this  is  not  known 
conclusively  at  present. 


9.2  One-dimensional  propagation 

The  assumption  of  one- dimensional  propagation  in  long  wave  theory  is 
quite  restrictive.  In  open  coastal  applications,  it  implies  a  complete  lack 
of  directionality  in  the  incident  wave  field.  The  results  of  various  studies 
(to  be  discussed  below)  have  indicated  that  results  of  one- dimensional 


model  calculations,  while  not  capable  of  reproducing  the  two-dimensional 
spatial  geometry  of  the  water  surface,  can  give  statistical  measures  which 
are  consistent  with  point  measurements  to  a  remarkable  degree.  For  the 
case  of  laterally  unbounded  motion  (or  for  motion  in  a  channel  of  constant 
width  26),  Johnson  [116]  derived  a  variable  depth  KdV  equation  which 
may  be  written  in  stationary  coordinates  as  [117] 

c  3c  ch? 

T]t  +  cqx  +  S-^T)  +  6—T]T]x  +  fi2—T)xxx  =  0(62,6fi2)  (141) 

Z  Zti  b 

where  h(x )  is  an  arbitary  function  of  x  with  a  derivative  of  0(8)  in 
size,  and  where  c  =  (6/6o)x^2  is  the  nondispersive  phase  speed  scaled 
by  a  value  based  on  the  reference  depth  h0.  This  equation  may  also 
be  altered  into  the  form  of  a  RLW  equation  by  changing  the  dispersive 
term.  Svendsen  and  Hansen  [118]  sought  to  solve  (141)  for  the  case  of 
cnoidal  waves  propagating  over  a  gently  sloping  bottom  by  assuming 
the  bottom  slope  to  be  an  order  smaller  than  indicated  above,  after 
which  a  first  order  solution  to  the  equation  gives  the  usual  cnoidal  wave 
as  the  local  solution.  The  second  order  correction  then  produces  an 
an  antisymmetric  perturbation  about  the  wave  crest,  leading  to  a  wave 
form  which  is  pitched  slightly  forward  at  the  crest.  It  is  not  clear  how 
applicable  an  approach  such  as  this  would  be  to  the  problem  of  arbitary 
topography,  since  several  phenomena  of  clear  importance  (such  as  the 
continued  deformation  of  a  shoaled  wave  over  a  shallow  shelf)  would  not 
be  reproduced.  Most  subsequent  attention  has  been  focussed  on  purely 
computational  approaches. 

For  waves  propagating  in  laterally-bounded  channels,  motion  becomes 
essentially  one-dimensional  if  the  horizontal  lengthscale  of  the  motion 
along  the  channel  greatly  exceeds  the  lengthscale  characterizing  the  chan¬ 
nel  width.  Peregrine  [119]  gave  a  method  for  constructing  Boussinesq 
equations  for  the  case  of  a  uniform  channel  of  arbitrary  cross-section,  and 
transformed  the  resulting  equations  to  a  KdV  equation  as  well.  Teng  h 
Wu  [120]  [121]  have  provided  extensions  of  both  a  Boussinesq  model  and 
a  KdV  model  to  the  case  of  a  channel  of  arbitrary  cross  section  which  is 
allowed  to  vary  along  the  channel  length.  For  the  case  of  a  rectangular 
channel  of  depth  h(x)  and  half  width  6(x),  the  KdV  form  of  their  model 
was  given  originally  by  Shuto  [122]  and  may  be  written  in  the  form 

,  ,  ^(bc)x  rZbc  , bch ?  „ 

brjt  +  bcrjx  +  S^-rj  +  b-—TjT]x  +  fi 2— —  rjxxx  =  0  (142) 
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Miles  [123]  has  pointed  out  that  (142)  is  non-conservative  of  mass,  with 
the  mass  conservation  equation  given  by 
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The  first  set  of  terms  in  brackets  in  (144)  represents  the  flux  of  mass 
across  a  station  x  —  xQ.  Miles  concludes  that  the  integral  term  in  (144) 
represents  a  transfer  of  mass  to  a  reflected  wave.  Kirby  &  Vengayil  [106] 
constructed  coupled  KdV  equations  for  opposite  going  waves  using  a 
splitting  technique,  verified  Miles’  conclusion,  and  further  established 
the  two-way  mechanism  for  mass  transfer  between  the  opposite  going 
waves.  Based  on  direct  computations  and  comparison  to  experimental 
data  [124],  Kirby  &  Vengayil  also  concluded  that  the  small  additional 
terms  that  Miles  used  to  rearrange  the  variable  coefficient  equation  into 
its  integral  form  should  be  include  in  the  model  system  from  the  outset, 
giving  a  rectangular  channel  equation  of  the  form 
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where  +(— )  denotes  rightward(leftward)  propagating  waves.  Mattioli  [125] 
has  further  extended  these  equations  by  including  the  nonlinear  coupling 
between  the  right  and  left  running  waves. 


10  Frequency  domain  formulations 

The  computation  of  a  time  dependent  wave  field  as  a  means  of  simulating 
a  complex  field  condition  presents  a  numerical  modeler  with  essentially 
the  same  problem  faced  by  the  experimentalist  -  how  to  interpret  the 
physical  implications  of  a  flood  of  data.  The  most  prevalent  data  analy¬ 
sis  technique  for  either  regular  waves  or  statistically  stationary  random 
waves  is  Fourier  analysis.  It  therefore  makes  some  sense  to  consider  a 
computation  of  the  desired  wave  field  in  the  frequency  domain  rather 
than  the  time  domain,  since  most  of  the  effects  being  looked  for  in  the 


frequency  domain  are  directly  available  from  the  computations  without 
further  analysis.  This  is  especially  true  for  shallow  water  waves,  where 
the  final  stage  of  highly  nonlinear  wave  evolution  can  be  characterized 
by  a  rapid  transfer  of  energy  from  the  fundamental  (or  spectral  peak 
in  the  case  of  random  waves)  to  higher  or  lower  harmonics,  by  means  of 
nearly  resonant  three- wave  interactions.  Much  of  the  original  impetus  for 
the  development  of  these  methods  came  from  a  study  of  the  water  wave 
analogy  of  the  process  of  second  harmonic  generation,  long  known  in  non¬ 
linear  optics  [126].  Mei  &  Unliiata  [127]  and  Boczar-Karakiewicz  [128] 
studied  the  phenomenon  of  second  harmonic  generation  and  recurrence 
in  weakly-dispersive  long  waves  using  the  spectral  approach,  and  showed 
that  the  method  provided  an  accurate  description  of  the  phenomenon; 
see  also  [129]  [130] 

In  this  section,  we  review  recent  results  obtained  using  spectral  meth¬ 
ods  in  the  context  of  the  variable  depth  problem. 


10.1  Deterministic  models  in  the  frequency  do¬ 
main 

Spectral  equations  for  wave  shoaling  and  evolution  have  been  presented 
in  a  number  of  different  formats  and  notations.  In  this  section,  we  will 
present  evolution  equations  for  complex  Fourier  amplitude  rather  than 
separating  the  resulting  equations  into  equations  for  real  amplitude  and 
phase. 

In  order  to  illustrate  the  structure  of  a  deterministic  model  for  water 
wave  shoaling,  we  derive  such  a  model  from  the  variable  depth  KdV 
equation  (141).  We  write  the  surface  displacement  rj  as  a  Fourier  series 
with  amplitudes  An  which  are  assumed  to  vary  on  a  slow  spatial  scale 
X  —  Sx  of  the  same  order  as  variations  in  water  depth; 

V  =  E  (146) 

where  ()*  denotes  a  complex  conjugation,  and  where  the  series  is  termi¬ 
nated  at  N  components  for  computational  reasons.  The  phase  function 
4>n  is  given  by 

V>(n)(x,<)  =  \  [  kn(X)dX  -  a >nt 
o  Jx 
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and  we  assume  that  the  wavenumber  and  frequency  for  each  model  are 
given  by  linear  nondispersive  theory  at  leading  order;  u2  =  k2h.  We  pro¬ 
ceed  by  substituting  (146)  in  (141),  collecting  terms  of  0(8,  fi2),  and  then 
isolating  the  contributions  to  the  motion  at  each  frequency  component 
n.  In  this  last  step,  we  assume  that  the  dominant  mechanism  leading  to 
nonlinear  interaction  between  different  frequency  components  is  resonant 
three  wave  interaction  (see  Phillips  [131],  Section  3.8).  We  thus  ignore 
any  contributions  from  quadratic  terms  which  are  not  oscillating  at  the 
same  frequency  as  the  target  frequency  component  n  and  obtain  the  set 
of  evolution  equations 


(148) 

The  resulting  model  equation  gives  the  rate  of  change  of  the  Fourier  am¬ 
plitude  for  each  frequency  component  due  to  the  leading-order  effects  of 
shoaling,  frequency  dispersion,  and  nonlinear  interaction.  This  model 
was  originally  given  by  Vengayil  &  Kirby  [132],  who  extended  it  to  in¬ 
clude  laminar  bottom  friction  effects  and  compared  it  to  data  for  shoaling 
periodic  waves  [63].  It  is  also  essentially  equivalent  to  the  “consistant” 
model  given  by  Freilich  &;  Guza  [133],  and  is  so  named  because  it  retains 
only  the  leading  order  effect  of  each  small  term  in  the  original  model 
equation.  Kirby  [134]  has  shown  that,  for  this  level  of  truncation,  the 
evolution  equation  being  solved  would  actually  be  most  analogous  to  a 
modified  KdV  equation  of  the  form 

3  1 

Vt  +  Vx-  6-r)T)t  -  =  0  (149) 

which  has  higher  wave  phase  speeds  and  broader  solitary  wave  crests 
than  corresponding  solutions  of  the  original  KdV  equations. 

Systems  of  the  form  (148)  may  be  applied  to  either  regular  or  irreg¬ 
ular  wave.  For  the  case  of  regular  wave  data,  it  is  sufficient  to  compute 
the  complex  Fourier  spectrum  An  for  one  period  of  the  wave,  and  use 
this  as  input  data  to  initialize  the  model  at  some  station  x.  The  equa¬ 
tions  are  then  solved  by  marching  in  x  to  obtain  the  complex  Fourier 
spectrum  at  more  shoreward  locations,  after  which  the  time  series  may 
be  reconstructed,  if  desired,  using  the  inverse  transform.  For  the  case  of 
irregular  waves,  the  process  used  is  essentially  similar  to  the  computation 
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of  smoothed  power  spectra  from  measurements.  The  time  series  of  data  is 
broken  into  a  set  of  shorter  series,  each  of  which  is  transformed  to  give  an 
ith  realization  of  complex  spectrum,  Aln.  The  individual  realizations  are 
then  used  to  initialize  separate  model  runs.  The  results  of  computations 
at  desired  x  locations  are  then  collected  and  used  to  compute  smoothed 
statistical  estimates.  For  example,  an  estimate  of  the  power  spectrum 
Sn  =  |An|2/2  based  on  /  realizations  of  the  data  would  be  given  by 

'  s.  =  ^EKI2  (iso) 
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This  process  can  be  somewhat  painstaking,  as  the  integration  of  (148) 
can  become  quite  stiff  and  cause  adaptive  step  size  schemes  to  grind  to 
a  halt.  Nevertheless,  it  is  a  reliable  approach  in  most  instances. 

Freilich  &  Guza  have  shown  that  a  model  equivalent  to  (148)  is  capa¬ 
ble  of  reproducing  most  of  the  features  of  spectral  evolution  during  the 
random  wave  shoaling  process,  including  amplification  of  the  spectrum 
do  to  shoaling  as  well  as  the  transfer  of  energy  to  higher  harmonics  of  the 
spectral  peak.  The  consistent  model  given  here  actually  tends  to  overpre¬ 
dict  purely  linear  shoaling  effects  due  to  the  lowest-order  representation 
of  the  shoaling  process  used.  The  shoaling  term  in  (148)  dictates  that 
the  amplitude  of  any  one  Fourier  component  evolves  (in  the  absence  of 
nonlinear  effects)  according  to  Green’s  Law 
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which  can  cause  serious  overprediction  of  amplitude  for  components  ini¬ 
tialized  at  a  sizeable  depth.  This  problem  is  partially  alleviated  by  re¬ 
taining  inhomogeneous  domain  effects  arising  from  the  dispersive  terms. 
Freilich  &  Guza  presented  a  model  which  retained  additional  effects  in 
both  shoaling  terms  and  nonlinear  coupling  coefficients;  this  model  has 
been  extensively  tested  against  data  and  has  been  seen  to  be  a  robust 
predictor  of  second  and  third-moment  statistics  [135]  [139],  even  in  cases 
where  wave  directional  spreading  becomes  significant. 

Madsen  &  Sprensen  [130]  have  considered  the  modification  of  the 
Boussinesq  model  for  the  purpose  of  obtaining  accurate  shoaling  coeffi¬ 
cients  in  spectral  models;  see  also  [45].  However,  since  each  frequency 
component  is  separately  identified,  it  is  in  fact  possible  to  represent 


full  linear  dispersive  effects  in  the  individual  equations.  Such  an  ap¬ 
proach  has  been  pursued  by  Agnon  et  al  [140],  who  used  the  spectral 
Zhakharov  equations  to  formulate  a  one-dimensional  model,  and  by  Kai- 
hatu  &  Kirby  [141],  who  developed  a  set  of  mild  slope  equations  coupled 
at  second-order  in  wave  steepness  and  extended  Agnon  et  al ’s  model  to 
two  horizontal  dimensions. 

For  the  case  where  directional  spreading  of  waves  becomes  significant, 
it  is  necessary  to  move  to  a  two-dimensional  or  weakly  two-dimensional 
format.  Liu  et  al  [103]  developed  a  two-dimensional  spectral  model  in  the 
form  of  a  coupled  set  of  parabolic  models;  this  model  is  reviewed  in  [68]. 
In  order  to  overcome  the  angular  restrictions  implied  by  the  use  of  a 
parabolic  approximation,  Kirby  [142]  developed  a  model  using  a  Fourier 
decomposition  in  both  time  and  the  longshore  direction  y.  The  model 
then  steps  a  complex  directional  spectrum  in  the  +a:-direction.  Several 
computational  examples  involving  laboratory  data  are  available  for  either 
of  these  models,  but  extensive  testing  against  field  data  is  lacking. 


10.2  Bispectra  and  third  moment  statistics 

As  waves  become  significantly  nonlinear  during  the  shoaling  process,  they 
become  asymmetric  in  both  the  vertical  and  horizontal  directions.  In 
terms  of  statistics,  the  vertical  asymmetry  is  referred  to  as  skewness 
and  is  given  by  the  normalized  third  moment  of  the  sea  surface,  while 
horizontal  asymmetry  is  called  asymmetry  and  is  given  most  directly  by 
the  normalized  third  moment  of  the  Hilbert  transform  of  the  sea  surface. 
Both  of  these  quantities  are  directly  related  to  a  third-order  spectrum 
known  as  the  bispectrum.  A  raw  complex  bispectral  value  is  given  in 
terms  of  complex  Fourier  amplitudes  by 

Bi,i  =  AiAjA’+j  (152) 

The  bispectrum  plays  a  central  role  in  the  detection  of  nonlinearity  in  a 
random  wave  train,  with  the  real  part  of  the  bispectrum  related  to  the 
presence  of  in-phase  harmonics  (as  in  a  regular  periodic  Stokes  or  cnoidal 
wave)  and  the  imaginary  part  related  to  the  presence  of  phase-locked  but 
out-of-phase  harmonics  (as  in  a  sawtooth  waveform).  In  this  guise,  the 
real  part  presents  the  the  distribution  of  contributions  to  skewness  over 
the  frequency  domain,  while  the  imaginary  part  gives  the  distribution  of 
contributions  to  asymmetry. 


In  the  context  of  spectral  models,  the  bispectrum  has  an  immediate 
dynamical  consequence.  If  we  rearrange  (148)  into  the  form  of  an  energy 
equation  (by  multiplying  by  A*  and  adding  An  times  the  conjugate  of 
(148)  ),  we  get  the  set  of  model  equations 
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where  5  denotes  the  imaginary  part.  The  presence  of  imaginary  bispec- 
tral  components  is  thus  directly  tied  to  the  transfer  of  energy  between 
spectral  components.  The  characteristic  pitched-forward  sawtooth  form 
of  surfzone  waves  is  representative  of  a  rapid,  sustained  transfer  of  energy 
towards  higher  frequencies.  This  transfer  must  continue  throughout  the 
surf  zone  as  long  as  the  waves  remain  bore-like  in  nature.  The  cessa¬ 
tion  of  breaking  and  reformation  of  a  less  asymmetric  individual  wave 
also  implies  a  cessation  of  this  energy  transfer  process.  This  result  has 
consequences  for  the  development  of  wave  breaking  models,  which  are 
considered  next. 


10.3  Wave  breaking 

The  problem  of  modelling  wave  breaking  in  a  spectral  model  is  fraught 
with  difficulties.  The  fact  that  dissipation  mechanisms  are  highly  local¬ 
ized  in  time  (i.e.,  confined  to  the  roller  region  as  the  roller  advects  past 
a  fixed  observer)  indicates  that  they  should  have  a  global  representation 
in  the  frequency  domain,  with  all  frequency  components  being  mutually 
dependent.  There  is  no  existing  model  which  includes  these  effects  at 
this  level  of  complexity.  Instead,  there  has  been  a  tendency  to  approach 
the  problem  by  appending  a  dissipation  term  to  the  spectral  model  (148), 
giving  a  model  of  the  form 

An,x  +  •  •  •  =  —ocnAn  (154) 

The  an  are  usually  chosen  by  using  a  bulk  energy  decay  model  to  compute 
the  energy  loss  for  the  entire  modelled  spectrum,  and  then  distributing 
this  loss  in  some  weighted  manner  across  the  frequency  spectrum.  An 
early  attempt  using  this  approach  was  made  by  Liu  [143],  who  studied 
the  breaking  of  regular  periodic  waves  [63].  Liu  chose  to  use  a  distribu¬ 
tion  of  ctn  s  with  no  dependence  on  frequency.  The  resulting  wave  forms 
in  the  modelled  surf  zone  show  a  marked  lack  of  asymmetry,  which  is  at 


odds  with  the  known  tendency  for  surf  zone  waves  to  become  sawtoothed 
in  form.  Based  on  the  arguments  in  the  previous  section,  it  is  clear  that 
this  approach  suppressed  the  transfer  of  energy  from  the  fundamental  to 
higher  harmonics,  indicating  an  over- accumulation  of  energy  at  relatively 
high  frequency.  Thus,  though  wave  height  and  total  variance  are  reason¬ 
ably  predicted,  third  moment  predictions  are  essentially  destroyed.  This 
approach  has  also  been  followed  in  a  recent  application  using  random 
waves  by  Eldeberky  &  Battjes  [144]. 

Mase  &;  Kirby  [83]  considered  the  problem  of  the  dependence  of  an 
on  frequency.  Based  on  observations  of  data  for  the  breaking  of  a  wave 
train  characterized  by  a  smooth,  Pierson-Moskowitz  spectrum,  they  con¬ 
cluded  that  an  should  have  a  quadratic  dependence  on  frequency.  Using 
this  distribution,  they  were  able  to  obtain  accurate  predictions  of  the 
evolution  of  wave  heights,  power  spectra  and  higher  moments  for  the 
simple  spectral  shapes  considered. 

Despite  the  indications  of  partial  success  observed  so  far,  the  accu¬ 
rate  prediction  of  transformation  and  breaking  of  complex,  multi-peaked 
spectra  has  still  not  been  accomplished.  It  is  likely  that  a  successful 
model  will  need  to  retain  some  of  the  global  nature  of  the  dissipation 
process  in  the  frequency  domain,  as  mentioned  above.  Models  of  this 
type  are  under  development  but  are  undocumented  at  present. 


10.4  Stochastic  model  formulations 

The  recent  implementation  of  deep  ocean  or  shelf  scale  wind  wave  models 
has  depended  to  a  large  extent  on  the  use  of  action  flux  models  aimed  at 
directly  computing  wave  statistics  as  opposed  to  any  phase-retaining  de¬ 
terministic  portrait  of  the  individual  waves  [145].  It  is  tempting  to  apply 
this  approach  in  the  nearshore  wave  environment  as  well,  although  the 
underlying  assumption  of  random  phases,  which  is  central  to  the  devel¬ 
opment  of  statistical  closures  in  intermediate  depth  action  flux  models,  is 
hard  to  justify  in  shallow  water  applications,  where  the  presence  of  rapid 
spectral  energy  transfer  implies  strong  phase  coupling  between  spectral 
components. 

The  framework  for  developing  stochastic  closures  for  spectral  models 
is  known  from  several  fields.  Equation  (153)  can  be  averaged  over  a  num¬ 
ber  of  realizations  to  yield  a  model  for  the  smoothed  power  spectrum  Sn 
which  depends  on  the  smoothed  bispectrum.  It  is  then  necessary  to  con- 


struct  a  model  to  compute  the  evolution  of  the  bispectrum  based  on  val¬ 
ues  of  teh  trispectrum,  or  fourth  order  spectral  moments.  This  cascade 
extends  to  all  orders,  and  hence  a  closure  problem  becomes  apparent. 
Work  in  this  format  is  in  its  infancy  and  will  not  be  well  documented  for 
some  time  yet.  There  is  some  indication  that  the  equations  for  the  bis¬ 
pectra  can  be  closed  in  terms  of  products  of  power  spectral  components. 
It  is  not  clear  that  this  closure  will  yield  a  simpler  computational  envi¬ 
ronment  than  the  deterministic  model  described  above,  as  the  number  of 
bispectral  components  essentially  goes  like  N2  and  thus  the  number  of 
equations  needed  to  compute  a  broad  spectrum  will  rival  the  number  of 
solutions  for  multiple  realizations  needed  in  the  deterministic  framework. 

There  is  one  existing  example  in  the  literature  of  an  attempt  to  in¬ 
voke  a  simpler  closure  and  thereby  confine  the  number  of  equations  to  be 
solved  to  the  number  of  spectral  components.  This  is  given  by  Abreu  et 
al  [146],  who  used  a  closure  proposed  by  Newell  &  Aucoin  [147]  which  is 
limited  to  strictly  resonant  interactions.  The  model  formulation  is  quite 
restrictive  in  that  it  neglects  all  interactions  between  waves  travelling  in 
different  directions,  as  is  required  by  the  limitation  to  tuned  resonances. 
However,  it  is  well  known  that  there  are  directional  interactions  in  spec¬ 
tral  wave  fields  [148]  [149].  Directional  interactions  are  also  crucial  in  the 
formation  of  hexagonal  patterns  in  biperiodic  wave  structures  and  in  the 
formation  of  Mach  stems  [142] . 

Despite  these  negative  conclusions  about  the  first  attempted  low- 
order  closure  scheme,  it  will  continue  to  be  highly  desireable  to  develop 
approximate  closures  which  eliminate  part  of  the  burden  of  computing 
the  entire  bispectrum  in  future  stochastic  models. 


11  Future  Directions 

The  maturation  of  the  theoretical  basis  and  numerical  methods  for  the 
Boussinesq  equations  has  led  to  a  recent  explosion  of  effort  aimed  at 
utilizing  these  models  for  coastal  engineering  prediction.  A  number  of 
purely  hydrodynamic  models  have  been  discussed  above.  Recently,  mod¬ 
els  coupling  the  hydrodynamics  to  sediment  transport  models  aimed  at 
predicting  seabed  evolution  have  also  been  developed  [150]  [151].  Further 
development  and  testing  in  this  area  needs  to  be  done. 

The  future  is  certain  to  see  an  extension  of  the  Boussinesq  and  similar 


equations  to  0(p4)  [152].  Since  the  present  0(p2)  models  with  corrected 
dispersive  properties  are  capable  of  predicting  the  propagation  of  waves 
over  most  of  the  range  of  shoaling  water  depths,  modifications  of  the 
model  equations  to  include  0(p4)  effects  are  likely  to  provide  a  limited  set 
of  improvements  if  propagation  characteristics  are  the  primary  concern. 
Indeed,  there  are  no  existing  data  sets  which  would  adequately  test  the 
improved  performance  of  an  0(p4)  propagation  model.  The  main  reason 
for  providing  such  an  extension  to  the  theory  would  be  to  provide  a 
more  accurate  predictive  capability  for  the  fluid  kinematics,  which  are 
presently  not  as  well  predicted  as  are  the  propagation  characteristics. 
Any  model  extension  in  this  direction  should  be  based  closely  on  a  related 
model  of  the  fluid  kinematics  carried  to  comparable  accuracy,  in  order  to 
reap  a  significant  benefit  from  the  extension. 

In  the  realm  of  spectral  applications,  the  future  is  likely  to  hold  the 
development  of  a  number  of  stochastic  models,  in  which  the  determin¬ 
istic  methods  prevalent  to  date  are  replaced  by  methods  for  computing 
the  wave  power  spectrum  together  with  higher  spectral  moments.  As 
mentioned  above,  these  methods  are  just  now  in  their  earliest  stages  of 
development. 
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